欠測値（欠損値）のを含むデータについて、データが完全にランダムに欠落しているかどうか、つまりMCAR（Missing
Completely at
Random）であるかを検定するLittleのMCAR検定をRで行う場合、やや古めのウェブページ等では{BaylorEdPsych}パッケージというすでにCRANからremoveされたパッケージのLittleMCAR()関数を使う方法を示しているものがある。
しかし、LittleのMCAR検定を行う場合、{BaylorEdPsych}パッケージのLittleMCAR()は基本的には使うべきではない*1。
これはCRANからremoveされたメンテナンスされてないパッケージ（および関数）だからというだけでなく、正しい統計量が計算できない場合が存在するからである。
LittleのMCAR検定を行うためには多変量正規分布の最尤推定を行う必要があるが、{BaylorEdPsych}パッケージのLittleMCAR()は{mvnmle}パッケージのmlest()関数を用いている。
このmlest()関数は規定値（mlest()の既定値というよりnlm()の既定値の関係で）では100回まで計算し、それまでに解が収束すればそれを出力し、100回まで解が収束しない場合はその時点での解を出力する仕様となっている。
そこで問題となるのは、mlest()は思いのほか既定値の範囲内では解が収束しない場合が多いということである。
例えばmlest()のドキュメントにあるexampleの例で扱われる13行5列のmissvalsデータですら、解が収束するまで311回を要する。
library(mvnmle) # iterationsが100となっている mlest(missvals)
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## $muhat
## [1] 5.279162 44.263596 10.349330 37.131345 94.477964
##
## $sigmahat
## [,1] [,2] [,3] [,4] [,5]
## [1,] 25.54034 13.547670 -23.684100 -13.87358 41.28384
## [2,] 13.54767 181.620955 -5.873532 -196.82020 144.87356
## [3,] -23.68410 -5.873532 55.856434 -39.06666 -61.74589
## [4,] -13.87358 -196.820204 -39.066663 272.76572 -119.67041
## [5,] 41.28384 144.873556 -61.745893 -119.67041 207.52799
##
## $value
## [1] 229.4776
##
## $gradient
## [1] -0.1398058 0.6927837 1.9804463 1.9536193 1.1419774 -5.1137371
## [7] 6.2012617 -9.6061015 -0.2843517 2.2540983 4.8672109 -4.3410681
## [13] -2.8165870 3.6930118 -1.6934243 5.6298682 3.8015192 7.0913520
## [19] 1.4127670 -2.1581259
##
## $stop.code
## [1] 4
##
## $iterations
## [1] 100
# 試しに（exampleの通り）iterlim = 400として計算回数の上限を上げてみる mlest(missvals, iterlim = 400)
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## $muhat
## [1] 6.655167 49.965258 11.769229 27.047091 95.423078
##
## $sigmahat
## [,1] [,2] [,3] [,4] [,5]
## [1,] 21.82557 20.86434 -24.900388 -11.473447 46.95304
## [2,] 20.86434 238.01242 -15.817386 -252.072282 195.60361
## [3,] -24.90039 -15.81739 37.869819 -9.599196 -47.55622
## [4,] -11.47345 -252.07228 -9.599196 294.182999 -190.59847
## [5,] 46.95304 195.60361 -47.556220 -190.598470 208.90486
##
## $value
## [1] 173.9566
##
## $gradient
## [1] 6.478684e-06 3.084000e-06 2.686993e-06 1.978164e-06 -1.767027e-06
## [6] -7.045778e-07 2.969584e-06 8.881344e-06 -4.205162e-05 1.562100e-05
## [11] -3.342251e-06 -1.177085e-05 7.482967e-05 6.159113e-06 2.863133e-05
## [16] -3.181941e-06 5.427961e-06 7.017551e-06 -4.230366e-06 -9.427307e-06
##
## $stop.code
## [1] 1
##
## $iterations
## [1] 311
もう一つ例を挙げると、組み込みデータで欠測値（欠損値）を含むairqualityもmlest()では同様に100回では収束しない。こちらは154回程度で収束する。
# iterationsが100になっている mlest(airquality)
## $muhat
## [1] 42.137122 185.930478 9.979428 77.870026 6.997105 15.801146
##
## $sigmahat
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1043.774622 898.167535 -65.2797197 209.522465 7.4784926 -6.8827520
## [2,] 898.167535 8052.761127 -17.0826526 232.911402 -7.3345985 -119.5057922
## [3,] -65.279720 -17.082653 12.3310937 -15.161552 -0.8830713 0.8089320
## [4,] 209.522465 232.911402 -15.1615518 89.016469 5.6073380 -10.8861366
## [5,] 7.478493 -7.334598 -0.8830713 5.607338 1.9932388 -0.1025824
## [6,] -6.882752 -119.505792 0.8089320 -10.886137 -0.1025824 78.0683925
##
## $value
## [1] 4641.691
##
## $gradient
## [1] -0.162113846 0.026192818 0.070248002 0.251939852 0.590048078
## [6] 0.050085350 -0.050411761 -0.077714705 0.018605852 -0.040893191
## [11] 0.036207894 -0.023638300 -0.003128662 -0.006800292 0.007657945
## [16] -0.048387847 -0.007636118 -0.570109478 0.001581611 -0.001259650
## [21] 0.068154804 -0.012631062 0.080040991 -0.009478754 1.124590199
## [26] 0.014077159 0.114649993
##
## $stop.code
## [1] 4
##
## $iterations
## [1] 100
mlest(airquality, iterlim = 200)
## $muhat
## [1] 42.167735 185.925880 9.980416 77.817784 6.989324 15.791587
##
## $sigmahat
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1043.974050 898.260489 -65.2747166 209.560887 7.46080107 -6.81740788
## [2,] 898.260489 8047.751854 -17.1140592 232.865524 -7.32823113 -119.43119444
## [3,] -65.274717 -17.114059 12.3341302 -15.178082 -0.88421910 0.84345576
## [4,] 209.560887 232.865524 -15.1780821 89.026436 5.60815570 -10.88566017
## [5,] 7.460801 -7.328231 -0.8842191 5.608156 1.99352220 -0.09948103
## [6,] -6.817408 -119.431194 0.8434558 -10.885660 -0.09948103 78.08865476
##
## $value
## [1] 4641.679
##
## $gradient
## [1] -1.009924e-01 2.628237e-02 3.929072e-03 -5.460913e-02 -5.338821e-03
## [6] -2.463106e-02 -8.983478e-02 1.252422e-01 -8.041472e-02 -2.128262e-02
## [11] -1.395802e-02 -9.496306e-02 -1.587978e-03 -2.612978e-03 1.055014e-04
## [16] 8.476491e-04 1.373337e-04 2.511115e-03 -3.738023e-04 -4.911271e-05
## [21] -1.971785e-03 9.822543e-04 1.492481e-03 -1.618901e-04 5.696165e-03
## [26] -3.336936e-03 4.702088e-04
##
## $stop.code
## [1] 2
##
## $iterations
## [1] 154
別に、{mvnmle}のmlest()上は正しく解が収束するまで計算するように指定してやればいい（ドキュメントのDetailsにnlm()つかってるからiterlim変えたらいいと書いてあるし、exampleも示している）だけだが、{BaylorEdPsych}パッケージのLittleMCAR()はデータフレーム（または行列）しか受け付けず、LittleMCAR()関数内で使用しているmlest()はデータフレームを与えるのみのため既定値である100回までしか計算しない設定になっているため、最尤推定するときに解が収束するまで計算していない結果を使用している可能性がある。
試しに{BaylorEdPsych}パッケージのLittleMCAR()関数の最尤推定の結果を取り出している箇所である、
gmean <- mlest(x)$muhat gcov <- mlest(x)$sigmahat
を次の通り、
gmean <- mlest(x, iterlim = 400)$muhat gcov <- mlest(x, iterlim = 400)$sigmahat
と計算回数の上限を400回程度に修正したものをLittleMCAR2()として、{BaylorEdPsych}パッケージのLittleMCAR()と結果を比較してみる。{BaylorEdPsych}はCRANからremoveされた当時の最新版の0.5を使用する。LittleMCAR()は結果にデータそのものを返すのが、それを掲載すると結果が長くなるので結果のうちの一部だけを載せる。
library(BaylorEdPsych) LittleMCAR(missvals)[1:4]
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## this could take a while
## $chi.square
## [1] 14.5923
##
## $df
## [1] 6
##
## $p.value
## [1] 0.02367611
##
## $missing.patterns
## [1] 3
LittleMCAR2(missvals)[1:4]
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## Warning in stats::nlm(lf, startvals, ...): NA/Inf
## は正の最大値で置き換えられました
## this could take a while
## $chi.square
## [1] 11.16069
##
## $df
## [1] 6
##
## $p.value
## [1] 0.08353528
##
## $missing.patterns
## [1] 3
二つの結果で統計量が変わっており、p値を見てみるとp = 0.05を基準として判断した場合、判断が変わってしまう。
もう一例、組み込みデータのairqualityの例を次に示す。
LittleMCAR(airquality)[1:4]
## this could take a while
## $chi.square
## [1] 35.14887
##
## $df
## [1] 14
##
## $p.value
## [1] 0.001397251
##
## $missing.patterns
## [1] 4
LittleMCAR2(airquality)[1:4]
## this could take a while
## $chi.square
## [1] 35.12481
##
## $df
## [1] 14
##
## $p.value
## [1] 0.001408772
##
## $missing.patterns
## [1] 4
この例ではp = 0.05を基準とした場合に判断そのものは変わらないが、やはり使った平均や共分散行列が異なっているので統計量の結果は異なっている。
このため、再度冒頭と同じことをいうがLittleのMCAR検定を行う場合、{BaylorEdPsych}パッケージのLittleMCAR()は基本的には使うべきではない。さらにいうと、CRANからremoveされメンテナンスされていないパッケージなので、この問題が今後解決する見込みは現在のところない。
解が収束していれば{BaylorEdPsych}パッケージのLittleMCAR()でも結果は同じだが、{BaylorEdPsych}パッケージのLittleMCAR()の結果で示される上記の情報からは、解が収束したのかどうか分からないというのも問題となる。
というわけで、現在としては、CRANに登録されている{naniar}パッケージのmcar_test()か、{misty}パッケージのna.test()を使ったほうが良い。ちなみに、{naniar}パッケージのmcar_test()も、{misty}パッケージのna.test()のどちらもドキュメントに、
https://web.archive.org/web/20201120030409/https://stats-bayes.com/post/2020/08/14/r-function-for-little-s-test-for-data-missing-completely-at-random/
というwebページを参考にしたとあり、このウェブページのコード自体は{BaylorEdPsych}パッケージのLittleMCAR()を参考にしている様子であり、基本的には同源である*2。
どちらも上記に示したLittleMCAR2()と結果が一致する（つまり解が収束された結果を用いている）。
library(naniar) mcar_test(missvals)
## # A tibble: 1 × 4
## statistic df p.value missing.patterns
## <dbl> <dbl> <dbl> <int>
## 1 11.2 6 0.0835 3
mcar_test(airquality)
## # A tibble: 1 × 4
## statistic df p.value missing.patterns
## <dbl> <dbl> <dbl> <int>
## 1 35.1 14 0.00142 4
library(misty)
## |-------------------------------------|
## | misty 0.5.4 (2023-11-14) |
## | Miscellaneous Functions T. Yanagida |
## |-------------------------------------|
na.test(missvals)
## Little's MCAR Test
##
## n nIncomp nPattern chi2 df pval
## 13 7 3 11.16 6 0.084
na.test(airquality)
## Little's MCAR Test
##
## n nIncomp nPattern chi2 df pval
## 153 42 4 35.11 14 0.001
こちらも、最尤推定の計算回数を指定する引数はないが、{naniar}パッケージのmcar_test()と{misty}パッケージのna.test()で最尤推定するために用いられている{norm}パッケージのem.norm()の計算回数の上限が既定値で1000回で多く、かつ{mvnmle}パッケージのmlest()よりも少ない回数で解が収束していることが多いため、さしあたって上記のmissvalsや、airquality程度のデータなら問題にならないためである。
library(norm)
## This package has some major limitations
## (for example, it does not work reliably when
## the number of variables exceeds 30),
## and has been superseded by the norm2 package.
s <- prelim.norm(as.matrix(missvals)) thetahat <- em.norm(s)
## Iterations of EM:
## 1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...
getparam.norm(s, thetahat)
## $mu
## [1] 6.655166 49.965259 11.769231 27.071654 95.423077
##
## $sigma
## [,1] [,2] [,3] [,4] [,5]
## [1,] 21.82557 20.86433 -24.900388 -11.597718 46.95303
## [2,] 20.86433 238.01244 -15.817377 -252.031412 195.60363
## [3,] -24.90039 -15.81738 37.869822 -9.400406 -47.55621
## [4,] -11.59772 -252.03141 -9.400406 293.791401 -190.77273
## [5,] 46.95303 195.60363 -47.556213 -190.772730 208.90485
途中に出力されるIteration of
EMというのが計算回数で124回で収束している（{mvnmle}パッケージのmlest()は311回）。
ただし、{norm}はパッケージを読み込んだときに出るメッセージの通り30以上の変数を含む場合に正しい結果を返さなない可能性があり、注意が必要である。
ところで、CRAN（のミラーサイトのうちPosit社が管理するもの）からパッケージがどの程度ダウンロードされているかを調べる{carnlogs}パッケージのcran_downloads()を使って、ここ1年間のダウンロード件数を調べると
library(cranlogs) cran_downloads("BaylorEdPsych", from = Sys.Date() - 365, to = Sys.Date()) |> summary()
## date count package
## Min. :2022-11-21 Min. : 0.000 Length:366
## 1st Qu.:2023-02-20 1st Qu.: 0.000 Class :character
## Median :2023-05-22 Median : 0.000 Mode :character
## Mean :2023-05-22 Mean : 1.016
## 3rd Qu.:2023-08-21 3rd Qu.: 1.000
## Max. :2023-11-21 Max. :41.000
となり、毎日1件くらいはダウンロードされ、日によっては41件もダウンロードされていることがわかり、世界の何処かでは今も使われているパッケージになっている（このうち1件は自分だけど）。
ここからは想像でしかないが、LittleのMCAR検定についてwebページで調べて、（代替案があるのにそんなはずはと思うかもしれないが）わざわざCRANのArchivesからダウンロードして使っている人も少なからずいるのではないかと思っている（他の関数を使っている例もあるとは思うが）。
ちなみにこの件は{mvnmle}パッケージのmlest()の挙動を眺めていて計算回数が既定値で解が収束しないパターンが多いことに気づいたことがきっかけで、{BaylorEdPsych}パッケージのLittleMCAR()のソースコードを読んで気づいた。意気揚々と書いてみたが、一箇所解けると割りと簡単に気づけそうな問題なので、他のところで書いている人がいるかも知れない。あるいは上記の内容で何か大きな勘違いをしているのであれば教えてほしい。
ついでにいうと、上記の内容に基づくと、過去に{BaylorEdPsych}パッケージ（{BaylorEdPsych}パッケージがCARNからremoveされる以前でも）のLittleMCAR()で計算した結果は統計量に誤りが合った可能性があると思う。
