欠損値（欠測値）を含むデータのmultivariate normal maximum-likelihood estimation、つまり多変量正規分布の平均ベクトルと共分散行列の最尤推定値の計算をRで行う方法。メモ。
CRANにあるパッケージでは{mvnmle}と{norm}で計算できる。ほかにも実現できるパッケージ等あるかもしれないが見つけられていない。
Jonathon D. Brown(2018). Advanced Statistics for the Behavioral Sciences A Computational Approach with R. Springer Cham. 354 - 356. [https://doi.org/10.1007/978-3-319-93549-2] で記載されているRでの実装版と比較する。
データは{mvnmle}パッケージの欠損値を含むappleデータを使う。
data("apple", package = "mvnmle")
{mvnmle}パッケージを使う
CRANから{mvnmle}をインストールして計算する。
# インストールしていない場合は # install.packages("mvnmle") # でインストールする library(mvnmle) mlest(apple)
## $muhat
## [1] 14.72227 49.33325
##
## $sigmahat
## [,1] [,2]
## [1,] 89.53415 -90.69653
## [2,] -90.69653 114.69470
##
## $value
## [1] 148.435
##
## $gradient
## [1] 4.988478e-06 2.892682e-06 8.726424e-07 1.682947e-05 -1.073488e-04
##
## $stop.code
## [1] 1
##
## $iterations
## [1] 34
ただし、ちょっと大きなデータを触るとすぐにnlm()がwornig
messageが出るが、計算結果は他の2例と変わらない。
# {mvnmle}のmissvalsデータをつかう data("missvals") # nlm()のiterlim（反復回数の上限回数）の既定値が100なので、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
# 次にとりあげる{norm}をつかって同じ計算をする s <- norm::prelim.norm(as.matrix(missvals)) thetahat <- norm::em.norm(s, showits = FALSE) norm::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
# 最後にとりあげるRでの実装での計算を行う。 emMLE(missvals$x1, missvals$x2, missvals$x3, missvals$x4, missvals$x5)
## $iter
## [1] 611
##
## $mu
## [1] 6.655166 49.965259 11.769231 27.047089 95.423077
##
## $sigma
## [,1] [,2] [,3] [,4] [,5]
## [1,] 21.82557 20.86433 -24.900388 -11.473441 46.95303
## [2,] 20.86433 238.01244 -15.817377 -252.072312 195.60363
## [3,] -24.90039 -15.81738 37.869822 -9.599211 -47.55621
## [4,] -11.47344 -252.07231 -9.599211 294.183044 -190.59849
## [5,] 46.95303 195.60363 -47.556213 -190.598491 208.90485
##
## $lnL
## [1] -41.0314
##
## $dif
## [1] 9.094947e-13
{norm}パッケージを使う
CRANから{norm}をインストールして計算する。
# インストールしていない場合は # install.packages("norm") # でインストールする 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(apple)) thetahat <- em.norm(s, showits = FALSE) getparam.norm(s, thetahat)
## $mu
## [1] 14.72222 49.33266
##
## $sigma
## [,1] [,2]
## [1,] 89.53395 -90.69081
## [2,] -90.69081 114.68298
ちなみに、{norm}関数の多くのドキュメントの参考文献が「See Section 5.3
of Schafer
(1994).」みたいないまいちわからない書き方がされているが、多分、mdataに参考文献として書いてあるSchafer,
J.L. (1997). Analysis of Incomplete Multivariate Data (1st ed.). Chapman
and Hall/CRC. https://doi.org/10.1201/9780367803025
が（目次しか見ていないが）本当の参考文献だと思う（関数ごとに年が微妙に異なるが）。
ついでにいうと、{norm}のmi.inference()の参考文献がp. 76-77 of Rubin
(1987)となっているが、これはRubin, D.B. (1987) Multiple Imputation for
Nonresponse in Surveys. John Wiley & Sons Inc., New York.
http://dx.doi.org/10.1002/9780470316696 のことだと思う。
さらにいうと、“This package … has been superseded by the norm2
package.”とパッケージを読み込んだときにメッセージが出るが、今年の8月に{norm2}はCRANからreomveされている。
（おまけ）{norm2}パッケージを使う
CRANからremoveされている{norm2}パッケージをつかって計算してみる。{remotes}パッケージを使ってインストールする。remotes::install_url("https://cran.r-project.org/src/contrib/Archive/norm2/norm2_2.0.4.tar.gz")とするとCRAN現在までに登録された最終版（removeされている）がインストールできる。
library(norm2) emNorm(apple) |> summary()
## Predictor (X) variables:
## Mean SD Observed Missing Pct.Missing
## CONST 1 0 18 0 0
##
## Response (Y) variables:
## Mean SD Observed Missing Pct.Missing
## size 14.72222 9.736563 18 0 0.00000
## worms 45.00000 10.539967 12 6 33.33333
##
## Missingness patterns for response (Y) variables
## (. denotes observed value, m denotes missing value)
## (variable names are displayed vertically)
## (rightmost column is the frequency):
## w
## so
## ir
## zm
## es
## .. 12
## .m 6
##
## Method: EM
## Prior: "uniform"
## Convergence criterion: 1e-05
## Iterations: 23
## Converged: TRUE
## Max. rel. difference: 9.3347e-06
## -2 Loglikelihood: 148.435
## -2 Log-posterior density: 148.435
## Worst fraction missing information: 0.6141
##
## Estimated coefficients (beta):
## size worms
## CONST 14.72222 49.33324
##
## Estimated covariance matrix (sigma):
## size worms
## size 89.53395 -90.69589
## worms -90.69589 114.69325
emNorm(missvals) |> summary()
## Predictor (X) variables:
## Mean SD Observed Missing Pct.Missing
## CONST 1 0 13 0 0
##
## Response (Y) variables:
## Mean SD Observed Missing Pct.Missing
## x1 6.00000 4.358899 9 4 30.76923
## x2 45.00000 15.953056 9 4 30.76923
## x3 11.76923 6.405126 13 0 0.00000
## x4 39.00000 16.492423 6 7 53.84615
## x5 95.42308 15.043723 13 0 0.00000
##
## Missingness patterns for response (Y) variables
## (. denotes observed value, m denotes missing value)
## (variable names are displayed vertically)
## (rightmost column is the frequency):
## xxxxx
## 12345
## ..... 6
## ...m. 3
## mm.m. 4
##
## Method: EM
## Prior: "uniform"
## Convergence criterion: 1e-05
## Iterations: 226
## Converged: TRUE
## Max. rel. difference: 9.5711e-06
## -2 Loglikelihood: 173.9567
## -2 Log-posterior density: 173.9567
## Worst fraction missing information: 0.9559
##
## Estimated coefficients (beta):
## x1 x2 x3 x4 x5
## CONST 6.655166 49.96526 11.76923 27.04733 95.42308
##
## Estimated covariance matrix (sigma):
## x1 x2 x3 x4 x5
## x1 21.82557 20.86433 -24.900388 -11.474685 46.95303
## x2 20.86433 238.01244 -15.817377 -252.071903 195.60363
## x3 -24.90039 -15.81738 37.869822 -9.597222 -47.55621
## x4 -11.47468 -252.07190 -9.597222 294.179111 -190.60023
## x5 46.95303 195.60363 -47.556213 -190.600234 208.90485
Rによる実装（Advanced Statistics for the Behavioral Sciences A Computational Approach with R.）
Rでの実装（Advanced Statistics for the Behavioral Sciences A
Computational Approach with
R.にある、emMLE()）と比較する。写経したコードを使うが、書籍のコードなので実装そのものはこの記事中はのせない。
なお、書籍はgoogle booksに登録されており、当該ページは読むことができる（現在のところ）。
emMLE(apple$size, apple$worms)
## $iter
## [1] 68
##
## $mu
## [1] 14.72222 49.33333
##
## $sigma
## [,1] [,2]
## [1,] 89.53395 -90.69673
## [2,] -90.69673 114.69496
##
## $lnL
## [1] -46.64932
##
## $dif
## [1] 6.110668e-13
計算速度の比較
appleデータを100回処理する時間を比較する。
# mvnmleの計算時間 mvnmle_time <- system.time(for(i in 1:100) mlest(apple)) # normの計算時間 norm_time <- system.time(for(i in 1:100){ s <- prelim.norm(as.matrix(apple)) thetahat <- em.norm(s, showits = FALSE) getparam.norm(s, thetahat) }) # Rでの実装の計算時間（emMLE） emMLE_time <- system.time(for(i in 1:100) emMLE(apple$size, apple$worms)) # norm2の計算時間 norm2_time <- system.time(for(i in 1:100) ans <- emNorm(apple)) # 計算結果をデータフレーム固める d <- list(mvnmle_time, norm_time, emMLE_time, norm2_time) |> Reduce(rbind, x = _) |> as.data.frame(row.names = NA) d <- cbind(data = c("mvnmle_time", "norm_time", "emMLE_time", "norm2_time"), d) d
## data user.self sys.self elapsed user.child sys.child
## 1 mvnmle_time 0.41 0.00 0.41 NA NA
## 2 norm_time 0.09 0.00 0.09 NA NA
## 3 emMLE_time 24.67 1.51 27.05 NA NA
## 4 norm2_time 0.08 0.00 0.08 NA NA
barplot(user.self ~ data, data = d)
計算の一部をCやFortranに投げている{mvnmle}、{norm}、{norm2}はRでの実装よりも早い。
少し大きなデータ（少し大きなデータと言ってもmissvalueはわずか13行5列のデータだが）を取り扱うとその差は更に大きくなる。
# mvnmleの計算時間 mvnmle_time <- system.time(for(i in 1:100) suppressWarnings(mlest(missvals, iterlim = 400))) # normの計算時間 norm_time <- system.time(for(i in 1:100){ s <- prelim.norm(as.matrix(missvals)) thetahat <- em.norm(s, showits = FALSE) getparam.norm(s, thetahat) }) # Rでの実装の計算時間（emMLE） emMLE_time <- system.time(for(i in 1:100) emMLE(missvals$x1, missvals$x2, missvals$x3, missvals$x4, missvals$x5)) # norm2の計算時間 norm2_time <- system.time(for(i in 1:100) ans <- emNorm(missvals)) # 計算結果をデータフレーム固める d <- list(mvnmle_time, norm_time, emMLE_time, norm2_time) |> Reduce(rbind, x = _) |> as.data.frame(row.names = NA) d <- cbind(data = c("mvnmle_time", "norm_time", "emMLE_time", "norm2_time"), d) d
## data user.self sys.self elapsed user.child sys.child
## 1 mvnmle_time 11.25 0.10 11.64 NA NA
## 2 norm_time 0.28 0.01 0.31 NA NA
## 3 emMLE_time 164.59 8.36 174.75 NA NA
## 4 norm2_time 0.34 0.00 0.35 NA NA
barplot(user.self ~ data, data = d)
{norm}と{norm2}がはやいが、{norm2}はCRANからはremoveされている。{norm}はパッケージを読み込んだときに出るメッセージの通り、潜在的な問題を抱えているらしい。
{mvnmle}は収束までの計算回数（iterations）が多い。あと、いろいろな条件を踏むと解がほかのパッケージで計算した結果やRで実装した方式で計算した結果とわずかに異なる場合があることに気づいた。もう少し調べていく。
Rでの実装は、小さなデータですら解が収束するまでの時間が結構かかる。またJonathon D. Brown(2018). Advanced Statistics for the Behavioral Sciences A Computational Approach with R. Springer Cham. 354 - 356. [https://doi.org/10.1007/978-3-319-93549-2] の実装だと、大きいデータだと計算できないことがある。
