井出草平の研究ノート

並列処理をするとBLRTのブートストラップは速くなるか

先日、潜在クラス分析のBLRT、要するにブートストラップを走らせるとやたらと時間がかかった。999回に設定してあることが大きな原因になっているが、Rががシングルコア処理をしているのも原因ではないかと考え、並列処理について検討してみた。MplusであればCPUのコア数を指定するコマンドがあるので、Rでもできないかと考え、試してみた。

ちなみに、先のエントリで実行したのは下記のスクリプトである。

library(randomLCA)
data(dentistry)
dentistry.lca2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=2)
dentistry.lca3 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=3)
obslrt <- 2*(logLik(dentistry.lca3)-logLik(dentistry.lca2))

## 以下BLRT
nsims <- 999
thesims <- simulate(dentistry.lca2, nsims)
simlrt <- rep(NA,nsims)
for (isim in 1:nsims) {
submodel <- refit(dentistry.lca2,newpatterns=thesims[[isim]])
fullmodel <- refit(dentistry.lca3,newpatterns=thesims[[isim]])
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel))
}

以下BLRTの所が時間がかかる場所。

処理時間の計測

処理時間は次のようにすると計測できるようだ。標準機能である。
stat.ethz.ch

ptm<-proc.time()

スクリプト

proc.time()-ptm

処理時間を知りたいスクリプトをサンドウィッチすると計測できる。

シングルコアでの計測

上記のBLRTの部分を走らせた結果が以下。。

ユーザ   システム       経過  
 272.84       0.75     273.61

ユーザーがRが実際に動いていたCPU時間で、経過が実際の経過時間なので273.61秒ということになる。

Microsoft R Open

mran.microsoft.com

Microsoft R OpenはCPUコア数を自動で判断して並列処理をする機能がある。

Microsoft R Openについてはこちらを参照のこと。

techtarget.itmedia.co.jp

RStudioでMicrosoft R Open 3.5.3を動かすと以下のようにCPU使用率が変わった。

f:id:iDES:20191122204603p:plain

いい感じで並列処理がされている。

処理時間は以下のようになった。

ユーザ   システム       経過  
1303.42       0.62     340.08

速くなっているのかなと思ったら、遅くなっていた。あれれ?

ユーザの時間の増加はマルチコア処理がされていてCPUの使用時間の総計が多くなっているためだろう。ユーザー時間が増えて、経過時間が遅くなっているということは、並列処理をしても早くならない、計算を仕分けしたり統合したりしていると、余計に時間がかかる、ことだろうか。

doSNOW

Rでは並列化させてた処理をする方法はいくつかある。代表的なdoSNOWパッケージを利用してみた。
下記のシミュレーションでもdoSNOWは良い値を出している。

itcweb.cc.affrc.go.jp

インストールして呼び出し、処理するスクリプトを下記のようにサンドイッチする。

atm<-proc.time()
library(doSNOW)
cl.SOCK <- makeCluster(12, type="SOCK")
registerDoSNOW(cl.SOCK)

スクリプト

stopCluster(cl.SOCK)
proc.time()-atm

結果は下記。

 ユーザ   システム       経過  
  290.06       1.09     295.94

結果は悪くないが、並列化をしないよりもやや遅い。

考察

結果失敗である。
スクリプトが並列化に対応すれば、処理も早くなるのだろうが、今の形では速くはならないようだ。
ベイズ統計学をするときにも並列化は必要になるが、RStanには並列化のコマンドがデフォルトで存在しているのでMicrosoft R OpenやdoSNOWは不要である。

ブートストラップは他の分析でも使用するが、今のところ僕は潜在クラス分析の時くらいしか使わないし、人生でそれほど数を走らせるわけではないので、時間とコストを考えると何も考えずシングルコアに頑張ってもらうのでもいいのかもしれない。

poLCAでデータ加工をしてRandomLCAでBLRTの計算をする[R]

私たちが使い慣れているデータ(行にケース、列に変数という形式のもの)を使ってRでブートストラップ尤度比検定(Bootstrapped Likelihood Ratio Test)を計算する方法についてのエントリである。
poLCAには残念ながらBLRTの機能がないので、poLCAのpredcell関数でデータを加工してからRandomLCAパッケージでBLRTの計算をするという手順を踏む。

データ

poLCAのcarcinomaデータを利用する。 子宮頸部のがんのデータである。デモなのでデータの説明は省く。

library(poLCA)
data(carcinoma)
carcinoma[100:105,] #100~105行目を抽出

次のようなデータである。

    A B C D E F G
100 2 2 2 2 2 1 2
101 2 2 2 2 2 1 2
102 2 2 2 2 2 1 2
103 2 2 2 2 2 2 2
104 2 2 2 2 2 2 2
105 2 2 2 2 2 2 2

predcell関数でデータ形式を変換

データ形式を下記の手順で変換する。

data(carcinoma)
f <- cbind(A,B,C,D,E,F,G)~1
clc2 <- poLCA(f,carcinoma,nclass=2)
carcinoma2<- clc2$predcell # スタック形式carcinoma2とする
print(carcinoma2) 

このようにするとRandomLCAで扱える形式のデータに変換できる。

   A B C D F G observed expected
1  1 1 1 1 1 1       36   30.028
2  1 2 1 1 1 1       10   16.197
3  1 2 1 1 1 2        6    2.023
4  2 1 1 1 1 1        2    3.751
5  2 1 2 1 1 2        1    0.203
6  2 2 1 1 1 1        4    2.023
7  2 2 1 1 1 2        8    4.075
8  2 2 1 1 2 2        1    2.769
9  2 2 1 2 1 2        3    4.447
10 2 2 1 2 2 2        3    3.222
11 2 2 2 1 1 2       13   11.858
12 2 2 2 1 2 2        5    8.592
13 2 2 2 2 1 2       10   13.797
14 2 2 2 2 2 2       16    9.997

もっとスマートな方法は他にもありそうだが、とりあえずpredcell関数しか僕は発見できていないので、今回はこの方法を利用する。

randomLCAパッケージで潜在クラス分析

データcarcinoma2のA~Gはデータは1/2で構成されている。RandomLCAは0/1である必要があるため、memiscパッケージを利用して0/1にリコードする。

library(memisc)
carcinoma2$A <- recode(carcinoma2$A , 0 <- 1, 1 <- 2)
carcinoma2$B <- recode(carcinoma2$B , 0 <- 1, 1 <- 2) 
carcinoma2$C <- recode(carcinoma2$C , 0 <- 1, 1 <- 2) 
carcinoma2$D <- recode(carcinoma2$D , 0 <- 1, 1 <- 2) 
carcinoma2$E <- recode(carcinoma2$E , 0 <- 1, 1 <- 2) 
carcinoma2$F <- recode(carcinoma2$F , 0 <- 1, 1 <- 2) 
carcinoma2$G <- recode(carcinoma2$G , 0 <- 1, 1 <- 2) 

データタイプは整数型でも実数型でも走るようだ。

今回は2クラスモデルと3クラスモデルを比較する。
比較する2クラスモデルと3クラスモデルを作成しておく。

library(randomLCA)
carcinoma2.lca2 <- randomLCA(carcinoma2[,1:7],freq=carcinoma2$observed,nclass=2)
carcinoma2.lca3 <- randomLCA(carcinoma2[,1:7],freq=carcinoma2$observed,nclass=3)

BLRT

BLRTのスクリプトに関する細かい説明は先のエントリを参照のこと。
まず、対数尤度の差の二乗。

## 3クラス、2クラスモデルの対数尤度の差の二乗
obslrt <- 2*(logLik(carcinoma2.lca3)-logLik(carcinoma2.lca2)) 
print(obslrt)

結果は以下。

'log Lik.' 47.09539 (df=23)

BLRT。

nsims <- 999 # ブートストラップの回数を指定
thesims <- simulate(carcinoma2.lca2, nsims) # シミュレーションを作成
simlrt <- rep(NA,nsims) 
for (isim in 1:nsims) {
submodel <- refit(carcinoma2.lca2,newpatterns=thesims[[isim]]) # submodelは比較対象の2クラスモデル
fullmodel <- refit(carcinoma2.lca3,newpatterns=thesims[[isim]]) # fullmodelは3クラスモデル
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel)) # 2つのモデルの尤度の差を二乗
}

P値。

print((sum(simlrt>=obslrt)+1)/(nsims+1)) # P値の計算

結果は以下。

p-value=0.001

Mplusに比べて手順が煩雑なのが欠点であろうか。

RandomLCAパッケージを用いた潜在クラス分析とBLRT[R]

Rで潜在クラス分析をするパッケージとしてはpoLCAが有名だが、RandomLCAでも潜在クラス分析ができる。 poLCAではブートストラップ尤度比検定(Bootstrapped Likelihood Ratio Test: BLRT)ができないが、RandomLCAではできるらしい。

cran.r-project.org

RandomLCAの特徴

RandomLCAは1) スタンダードな潜在クラス分析 2)シングルレベルランダム効果モデル 3)2レベルランダム効果モデル3つが使用できる。

ちなみにこのパッケージの作者はBeath自身である。 RandomLCAは他の統計パッケージではできない解析ができるようだ。

ランダム効果モデルの話はまたの機会にするとして、今回はとりあえず、スタンダードな潜在クラス分析を行ってみたい。

データの特徴

RandomLCAで利用するデータの形式は私たちが普段使っている形式とは異なっている。サンプルデータのHIV testing dataを例にとろう。

f:id:iDES:20191119174243p:plain

V1~V4までの0/1のパターンは2の4乗なので16通りあり、右側にその条件を満たす頻度が書かれている。もともとは潜在クラス分析の計算はこの形式にしてから計算をするので、慣れない感じはあるが、こちらのデータセット形式の方が正統的ではある。このブログでは一度も言及したことはずだが、統計パッケージのLEMもこのような形式のデータセットを用いる。

潜在クラス分析

サンプルデータとして利用するのはHandelman et al. (1986)の歯のX線データである。 1-5行目が変数であり、6行目はfreqという頻度のデータである。

library(randomLCA)
data(dentistry) 
dentistry.lca2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=2)
summary(dentistry.lca2)

結果は次のようになる。

  Classes      AIC      BIC     AIC3    logLik penlogLik
        2 14952.77 15021.64 14963.77 -7465.385 -7465.446
Class probabilities 
Class  1 Class  2 
  0.1961   0.8039 
Outcome probabilities 
             V1     V2     V3     V4     V5
Class  1 0.4034 0.7129 0.5981 0.4888 0.9155
Class  2 0.0106 0.1020 0.0136 0.0316 0.3053

BLRT

BLRTは下記のように実施する。 今回は2クラスモデルと3クラスモデルを比較する。

dentistry.lca2 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=2) # 2クラスモデル
dentistry.lca3 <- randomLCA(dentistry[,1:5],freq=dentistry$freq,nclass=3) # 3クラスモデル
obslrt <- 2*(logLik(dentistry.lca3)-logLik(dentistry.lca2)) # 3クラス、2クラスモデルの対数尤度の差の二乗
print(obslrt)

Printはこの出力は、観察された(ブートストラップでシミュレーションされたものではない)実際の対数尤度の差の二乗である。この値はP値を求める際に使用する。

'log Lik.' 108.3152 (df=17)

BLRTは下記の次のように行う。

nsims <- 999 # ブートストラップの回数を指定
thesims <- simulate(dentistry.lca2, nsims) # シミュレーションを作成
simlrt <- rep(NA,nsims) 
for (isim in 1:nsims) {
submodel <- refit(dentistry.lca2,newpatterns=thesims[[isim]]) # submodelは比較対象の2クラスモデル
fullmodel <- refit(dentistry.lca3,newpatterns=thesims[[isim]]) # fullmodelは3クラスモデル
simlrt[isim] <- 2*(logLik(fullmodel)-logLik(submodel)) # 2つのモデルの尤度の差を二乗
}

P値はシミュレートされたLRTが実際のLRTより大きい割合である。1は修正値である。

print((sum(simlrt>=obslrt)+1)/(nsims+1)) # P値の計算

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

0.001

P値が有意であった場合、2クラスモデルから3クラスモデルにしてモデルが改善した、つまり2ではなく3クラスモデルが支持されるという解釈を行う。

課題と考察

P値の計算式の分母に反復回数であるnsimsが入っているので、P値を細かく出そうとすると反復回数はそこそこ大きい数字の方が求められる。小数点3桁でP値を表示するには999+1回の反復が必要である。

この例ではシミュレートされたLRTが実際のLRTより大きかった回数(sum(simlrt>=obslrt)+1)は1回もないので値は0である。Rでsum(simlrt>=obslrt)+1と打って返されるのは0であることで確認できる。そこに修正値の1を足しているので分子は1である。

反復回数が9回だった場合、このモデル比較だと、P値は0+1/(9+1)であり、0.1、つまり10%となる。頻度主義は5%以下をだいたいの場合は有意とするので、モデルの改善がないと判断される。つまり、BLRTでは反復回数は完全に自由に決めてよいというわけではないのだ。

修正値が1の場合は反復回数は少なくとも20回ないと5%有意にならない。20回というのは分子が0(クラス数の増加によつて明らかにモデル改善がある場合)だけであって、ぎりぎりの数値であり、実際には20回というのは少なすぎる。

修正値を0.1にすると、反復回数は100回で小数点3桁のP値が得られる。修正値は分子・分母(特に分母)を0にしないために入れられるものなので、1である必要はないのかもしれない。0.1であってもよさそうに思える。ブートストラップの反復回数と修正値は連動しているというのは少なくとも間違いないだろう。修正値として何が適切なのかを知るには、BLRTについてちゃんと勉強しないといけないようである。

潜在クラス分析のエントロピーの計算[R]

RのpoLCAパッケージで潜在クラス分析のエントロピーを計算してみたい。

相対的エントロピーと絶対的エントロピー

エントロピーという指標はMplusユーザーにとってはお馴染みのものである。
Mplusで出力されるのは正確に言うと、相対的エントロピー(Relative Entropy)と呼ばれるものである。RのpoLCAパッケージでもエントロピーが出力できるが、絶対的エントロピーであって、相対的エントロピーとは異なったものだ。

相対的エントロピーの文献はRamaswamyらの論文であろうか。

Ramaswamy et al. 1993, "An Empirical Pooling Approach for Estimating Marketing Mix Elasticities with PIMS Data"

相対化エントロピーは最小値が0で1が最大値である。1に近い方がクラスの分類が正確だとされている。
相対化エントロピーの値は0.7以上で十分、0.8以上で良好な分類だと言われている。
Ramaswamyらの論文ではこの基準が書かれているかというと、似たような記述はあるが、決して基準の提唱をしているわけではない。Mplusの作者Muthen & MuthenがMplusのマニュアルで0.8以上の推奨をしていることが確認できた。

下記の本の304頁には0.7以上で十分、0.8以上で良好という記述がある。

International Perspectives on Teacher Knowledge, Beliefs and Opportunities to Learn: TEDS-M Results (Advances in Mathematics Education)

International Perspectives on Teacher Knowledge, Beliefs and Opportunities to Learn: TEDS-M Results (Advances in Mathematics Education)

該当の論文はResearch Gateにある。

Sigrid Blömeke and Gabriele Kaiser - Homogeneity or Heterogeneity? Profiles of Opportunities to Learn in Primary Teacher Education and Their Relationship to Cultural Context and Outcomes

本来は異なったクラス数のモデル間で、相対化していないエントロピーを比較するのが正しい使い方なのではないかと思うが、標準化してもこの機能は果たせるので、相対化エントロピーでも問題はない。ただ、0.7以上といった基準は個人的にはあまりどうでもいい気がする。この議論は相関係数がどの程度あれば強い相関なのか、といった議論と同種のものであるからだ。

絶対的エントロピーの計算

潜在クラス分析は先のエントリーで使用したものを利用する。先のエントリでは共変量を設定しているが下記のスクリプトでは共変量のないノーマルな潜在クラス分析を行っている。

library("poLCA")
data("election",package="poLCA") # electionデータ読み込み
election$GENDER <- factor(election$GENDER) # 因子型にする
d1 <- election[,c(1:3,7:12,16:17)] # 道徳性・思いやり・知識が豊富のゴア、ブッシュの値、性別、政党を抜き出し
d1 <- na.omit(d1) # 欠損値除去
fm <- cbind(MORALG,CARESG,KNOWG,MORALB,CARESB,KNOWB)~1 # 式formulaの作成
res3 <- poLCA(formula=fm,data=d1,nclass=3,nrep=10) # 潜在クラス分析の実行
poLCA.entropy(res3) # エントロピーの計算

エントロピーは6.62と計算される。

6.621753

この関数でエラーが出るという報告がある(こちら)。僕の環境ではエラーは出ていないので、特に問題ないのだが、エラーが出る人もいるのだろうか。

その代わりに次のような手計算の式が書かれてあるのだが、僕の場合はこちらの計算でエラーが出た。

poLCA.entropy.fix <- function (lc)
{
 K.j <- sapply(lc$probs, ncol)
 fullcell <- expand.grid(sapply(K.j, seq, from = 1))
 P.c <- poLCA.predcell(lc, fullcell)
 return(-sum(P.c * log(P.c), na.rm = TRUE))
}
poLCA.entropy.fix(lca.polca)
[1] 3.602397

相対的エントロピーの計算

RのpoLCAは相対的エントロピーの計算はしてくれないので、自力で行う必要がある。自力といっても、潜在クラス分析の結果はもちろん使って行うので、先に潜在クラス分析の計算はしなくてはならない。下のスクリプトこちらを参考にした。

n <- 1346 # サンプルサイズ
K <- 3 #クラス数
nume.E <- -sum(res3$posterior * log(res3$posterior)) # 分子
deno.E <- n*log(K) #分母
Entro <- 1-(nume.E/deno.E)
Entro

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

0.7329062

他の計算法としては、Daniel Oberskiが提案しているものがある(PDF)。

entropy <- function (p) sum(-p*log(p))
error_prior <- entropy(res3$P) # クラス割合
error_post <- mean(apply(res3$posterior, 1, entropy))
R2_entropy <- (error_prior - error_post) / error_prior
R2_entropy

共変量を伴った潜在クラス分析[R]

Rで共変量を伴った潜在クラス分析を行ってみたい。
パッケージはpoLCAを使う。

poLCAパッケージを利用した潜在クラス分析は以前にエントリしている。
RのpoLCAパッケージで潜在クラス分析を行う

共変量とは何かというのは、以前にMplusでの方法をエントリしているのでそちらを参照のこと。
Mplus 共変量を伴った潜在クラス分析

パッケージを読み込む。

library("poLCA")

データはpoLCAパッケージに含まれるelectionデータを利用する。 2000年のアメリカ大統領選挙アル・ゴアジョージ・W・ブッシュの評価のデータである。

下記の6点について4段階で評価している。

  1. 道徳性 moral
  2. 思いやり caring
  3. 知識が豊富 knowledgeable
  4. 良いリーダー good leader
  5. 正直ではない dishonest、
  6. 知的 intelligent

他の変数としては以下のものが含まれる。

  • 投票 VOTE3:(1) ゴア (2) ブッシュ (3) その他
  • 年齢 AGE: 連続変数
  • 教育 EDUC (1)第8学年以下 (2)9-11年生(3)高校卒業または同等(4)12年以上の学校教育(5)ジュニアまたはコミュニティカレッジレベルの学位 (6)学士レベルの学位 (7)高度な学位
  • 性別 GENDER:(1)男性 (2)女性
  • 政党 PARTY:(1)強い民主党員(2)弱い民主党員 (3)所属の無い民主党員 (4)支持政党なし(5)所属の無い共和党員 (6)弱い共和党員 (7)強い共和党

VOTE3を用いれば、抹消結果Distal Outcomesを伴った潜在クラス分析もできるデータのようだ。

少しデータが多いので、分析するデータを編集。

data("election",package="poLCA") # electionデータ読み込み
election$GENDER <- factor(election$GENDER) # 因子型にする
## 道徳性・思いやり・知識が豊富のゴア、ブッシュの値、年齢、教育、性別、政党を抜き出し
d1 <- election[,c(1:3,7:12,14:17)]
d1 <- na.omit(d1) # 欠損値除去

共変量を伴った潜在クラス分析

共変量は式の最後のチルダの後につける。
ノーマルな潜在クラス分析だとここに~1と表記するが、そこに共変量を設定することになる。 クラス数は3で走らせる。

fm1 <- cbind(MORALG,CARESG,KNOWG,MORALB,CARESB,KNOWB)~PARTY # 式formulaの作成
res1 <- poLCA(formula=fm1,data=d1,nclass=3,nrep=10) # 潜在クラス分析の実行

僕の環境ではRStudioのR Markdownからは出力ができず、コンソールからのみ共変量の結果は出力できた。
それにしても、よく起こるバグである。

共変量のところのみ結果を表示する。

========================================================= 
Fit for 3 latent classes: 
========================================================= 
2 / 1 
            Coefficient  Std. error  t value  Pr(>|t|)
(Intercept)     2.44815     0.21781   11.240         0
PARTY          -0.74290     0.06139  -12.101         0
========================================================= 
3 / 1 
            Coefficient  Std. error  t value  Pr(>|t|)
(Intercept)    -3.88769     0.47323   -8.215         0
PARTY           0.66537     0.08135    8.179         0
========================================================= 

p値は0である。ゴアとブッシュへの評価が潜在クラスの中身なので、政党支持が関連を持っているのは明らかだろう。
参照カテゴリは1である。参照カテゴリを変更させるコマンドは探してみたがなかった。変更できないと少し不便だ。

2つ以上の共変量を伴った潜在クラス分析

2つ以上の共変量を設定するときは単純に+をいれるとよい。

fm2 <- cbind(MORALG,CARESG,KNOWG,MORALB,CARESB,KNOWB)~PARTY+AGE+EDUC # 式の作成
res2 <- poLCA(formula=fm2,data=d1,nclass=3,nrep=10) # 潜在クラス分析の実行

今回も共変量のところのみ結果を表示する。

========================================================= 
Fit for 3 latent classes: 
========================================================= 
2 / 1 
            Coefficient  Std. error  t value  Pr(>|t|)
(Intercept)     4.38411     0.65075    6.737     0.000
PARTY          -1.59906     0.10977  -14.567     0.000
AGE             0.00910     0.00825    1.102     0.270
EDUC            0.27591     0.08981    3.072     0.002
========================================================= 
3 / 1 
            Coefficient  Std. error  t value  Pr(>|t|)
(Intercept)     5.02711     0.58119    8.650     0.000
PARTY          -0.99110     0.09064  -10.935     0.000
AGE            -0.02268     0.00714   -3.179     0.002
EDUC            0.16649     0.07637    2.180     0.029
========================================================= 

潜在クラス分析の応答確率の棒グラフ[Stata]

Stataでは潜在クラス分析の応答確率の棒グラフが描ける。棒グラフというのは今まで見たことがなかった。不要な気もしないでもないが、説明するときにわかりやすいのかもしれない。

使用するのはmarginsmarginsplotである。説明によると条件付き応答確率の関数であるestat lcmeanmarginsのラッパーらしい。ちなみにestat lcmean前のエントリで応答確率を計算するときに使ったコマンドである。

前のエントリで計算はしているが、一応今回も同じ潜在クラス分析のモデルを作っておこう。

. use http://www.stata-press.com/data/r16/gsem_lca1.dta, clear
. gsem (accident play insurance stock <- ), logit lclass(C 2)

前回のフィッテング指標の結果2クラスの方が良かったので、今回は2クラスのモデルで計算してある。プロットは下記のように記述する。

. margins, predict(outcome(accident) class(1)) ///
predict(outcome(play) class(1)) ///
predict(outcome(insurance) class(1)) ///
predict(outcome(stock) class(1))

. marginsplot, recast(bar) title("Class 1") xtitle("") ///
xlabel(1 "accident" 2 "play" 3 "insurance" 4 "stock", angle(45)) ///
ytitle("Predicted mean") ylabel(0(0.1)1) name(class1)


. margins, predict(outcome(accident) class(2)) ///
predict(outcome(play) class(2)) ///
predict(outcome(insurance) class(2)) ///
predict(outcome(stock) class(2))

. marginsplot, recast(bar) title("Class 2") xtitle("") ///
xlabel(1 "accident" 2 "play" 3 "insurance" 4 "stock", angle(45)) ///
ytitle("") ylabel(0(0.1)1) name(class2)


.graph combine class1 class2, cols(2)

書いてあることは難しくないが、内容のわりにやたらと長い。
2クラスなのでこの程度だが、4クラスだとこの倍の量になる。とても洗練されているとは言い難い。estat lcmeanmarginsのラッパーであれば、プロットもオプションとして用意して、marginsから記述しなくてもいいようにできなかったのだろうか。

少し注意が必要そうなのは、ylabelの目盛りを調節だ。この図を出すまでに何回か調節をする必要があった。

f:id:iDES:20191116032605p:plain