セットアップ
本日使用するパッケージとLMSからダウンロードして実習用データを読み込む。読み込むデータはdfに格納し、中身を確認してみよう。
library (tidyverse)library (summarytools)library (marginaleffects)<- read_csv ("Data/M1.csv" )
# A tibble: 60 × 8
No Year Name Duration First Final Order Score_Mean
<dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 13 2017 ゆにばーす 5 1 0 1 89.4
2 13 2017 カミナリ 7 0 0 2 88.3
3 13 2017 とろサーモン 16 1 1 3 92.1
4 13 2017 スーパーマラドーナ 15 1 0 4 91.4
5 13 2017 かまいたち 14 0 0 5 91.4
6 13 2017 マヂカルラブリー 11 1 0 6 86.7
7 13 2017 さや香 4 1 0 7 89.7
8 13 2017 ミキ 6 1 1 8 92.9
9 13 2017 和牛 12 0 1 9 93.3
10 13 2017 ジャルジャル 15 0 0 10 90.9
# … with 50 more rows
続いて、dfの記述統計量を計算する。コンビー名（Name）列の記述統計は不要であるため、予め除外しておこう。
%>% select (- Name) %>% descr (stats = c ("mean" , "sd" , "min" , "max" , "n.valid" ),transpose = TRUE , order = "p" )
Descriptive Statistics
df
N: 60
Mean Std.Dev Min Max N.Valid
---------------- --------- --------- --------- --------- ---------
No 15.50 1.72 13.00 18.00 60.00
Year 2019.50 1.72 2017.00 2022.00 60.00
Duration 10.60 3.70 2.00 16.00 60.00
First 0.52 0.50 0.00 1.00 60.00
Final 0.30 0.46 0.00 1.00 60.00
Order 5.50 2.90 1.00 10.00 60.00
Score_Mean 91.33 2.40 84.86 97.29 60.00
データの紹介は以下ととおりである。
No第X回大会
Year大会年度
Nameコンビ名
Duration結成からの経過年数
First初出場ダミー
Finalファイナルステージへの進出有無
Order出場順番
1から10
Score_Mean平均得点
7人の審査委員からの評価の平均値
今回はdfを使って、M-1グランプリにおける出場順番と得点の関係を調べるとする。応答変数は得点の合計を審査委員の数で割った値である平均得点（Score_Mean）とファイナルステージへの進出有無（Final）である。説明変数は出場順番（Order; 主な説明変数）、初出場ダミー（First）、コンビー結成からの経過年数（Duration）である。また、出場順番と初出場ダミー、または経過年数の交互作用も検討し、交互作用を仮定しない変数は統制変数として投入する。ここでは以下の4モデルを推定するとしよう。
fit_lm: 線形回帰分析（応答変数: 平均得点）fit_logit: ロジスティック回帰分析（応答変数: ファイナルステージへの進出有無）fit_inter1: 線形回帰分析 + 交互作用（ダミー変数）
fit_inter2: 線形回帰分析 + 交互作用（連続変数）
<- lm (Score_Mean ~ Order + First + Duration, data = df)<- glm (Final ~ Order + First + Duration, data = df, family = binomial ("logit" ))<- lm (Score_Mean ~ Order * First + Duration, data = df)<- lm (Score_Mean ~ Order * Duration + First, data = df)
{modelsummary}のmodelsummary()関数を使って、4つのモデルを概観してみよう。
modelsummary (list ("線形" = fit_lm, "ロジスティック" = fit_logit, "交互作用（1）" = fit_inter1, "交互作用（2）" = fit_inter2),statistic = "({p.value})" )
線形
ロジスティック
交互作用（1）
交互作用（2）
(Intercept)
89.605
−3.121
89.825
88.501
(<0.001)
(0.036)
(<0.001)
(<0.001)
Order
0.258
0.276
0.210
0.461
(0.016)
(0.025)
(0.189)
(0.192)
First
−0.963
−1.126
−1.421
−1.004
(0.128)
(0.092)
(0.276)
(0.117)
Duration
0.076
0.109
0.078
0.176
(0.381)
(0.251)
(0.370)
(0.348)
Order × First
0.085
(0.687)
Order × Duration
−0.019
(0.544)
Num.Obs.
60
60
60
60
R2
0.144
0.146
0.150
R2 Adj.
0.098
0.084
0.088
AIC
275.0
71.2
276.8
276.6
BIC
285.5
79.6
289.4
289.2
Log.Lik.
−132.495
−31.618
−132.405
−132.292
F
3.134
2.547
2.356
2.417
RMSE
2.20
0.42
2.20
2.19
線形回帰分析
線形回帰分析の場合、推定結果を表で示すだけでも良いが、予測値のグラフも示すと読者によってより理解しやすくなる。今回の例（fit_lm）の場合、出場順番（Order）と平均得点間の関係をグラフを示してみよう。出場順番は1番目から10番目まであるため、Orderの値が1、2、3、…、10の場合に対し、予測値（平均得点の予測値）とその95%信頼区間を出してみよう。使用するのは{marginaleffects}のpredictions()関数である。
第1引数は回帰分析のオブジェクト名、第2引数はnewdataにdatagrid()関数を割り当てたものである。datagrid()の中には説明変数名 = 値を指定する。1のみならOrder = 1の場合の予測値が計算され、c(1, 2, 3)ならOrderの値が1、2、3の場合の予測値が計算される。公差1の等差数列は:演算子で表現することもできるので、今回は1:10とする。計算結果はlm_predに格納し、出力する。
<- predictions (fit_lm, newdata = datagrid (Order = 1 : 10 ))
rowid type predicted std.error statistic p.value conf.low conf.high
1 1 response 90.16668 0.5527105 163.1355 0 89.05947 91.27389
2 2 response 90.42488 0.4679791 193.2242 0 89.48740 91.36235
3 3 response 90.68308 0.3926080 230.9761 0 89.89659 91.46957
4 4 response 90.94127 0.3330148 273.0848 0 90.27417 91.60838
5 5 response 91.19947 0.2987938 305.2254 0 90.60092 91.79803
6 6 response 91.45767 0.2987938 306.0896 0 90.85911 92.05623
7 7 response 91.71587 0.3330148 275.4108 0 91.04876 92.38298
8 8 response 91.97407 0.3926080 234.2643 0 91.18758 92.76055
9 9 response 92.23226 0.4679791 197.0863 0 91.29479 93.16974
10 10 response 92.49046 0.5527105 167.3398 0 91.38325 93.59767
Score_Mean First Duration Order
1 91.32857 0.5166667 10.6 1
2 91.32857 0.5166667 10.6 2
3 91.32857 0.5166667 10.6 3
4 91.32857 0.5166667 10.6 4
5 91.32857 0.5166667 10.6 5
6 91.32857 0.5166667 10.6 6
7 91.32857 0.5166667 10.6 7
8 91.32857 0.5166667 10.6 8
9 91.32857 0.5166667 10.6 9
10 91.32857 0.5166667 10.6 10
大量の列が出力されているが、作図のために必要な変数はOrder（出場順番）、predicted（予測値）、conf.low（95%信頼区間の下限）、conf.high（95%信頼区間の上限）列のみだ。これらにstd.error列を加え、lm_predを上書きしよう。
<- lm_pred %>% select (Order, predicted, std.error, conf.low, conf.high)
Order predicted std.error conf.low conf.high
1 1 90.16668 0.5527105 89.05947 91.27389
2 2 90.42488 0.4679791 89.48740 91.36235
3 3 90.68308 0.3926080 89.89659 91.46957
4 4 90.94127 0.3330148 90.27417 91.60838
5 5 91.19947 0.2987938 90.60092 91.79803
6 6 91.45767 0.2987938 90.85911 92.05623
7 7 91.71587 0.3330148 91.04876 92.38298
8 8 91.97407 0.3926080 91.18758 92.76055
9 9 92.23226 0.4679791 91.29479 93.16974
10 10 92.49046 0.5527105 91.38325 93.59767
だいぶ読みやすい表としてまとまった。あとは、このlm_predオブジェクトを使って作図をするだけだ。点推定値と区間を同時に示すときにはpoint-rangeプロットが有効である（本資料の後半では線と面を組み合わせたグラフの作り方も紹介するが、その方法でも良い）。使用する幾何オブジェクトはgeom_pointrange()であり、aes()内部でマッピングを行う。最低でもx、y、ymin、ymaxのマッピングが必要であり、xとyはそれぞれ横軸と縦軸に対応する変数、yminとymaxは区間の下限と上限に対応する変数を指定する。また、横軸の目盛りを1刻みにするため、scale_x_continuous()レイターを足す。scale_x_continuous()はbreaksに目盛りの位置を、labelsには各目盛りに付けるラベルのベクトルを割り当てる。breaksとlabelsは同じ長さのベクトルである必要がある。
%>% ggplot () + geom_pointrange (aes (x = Order, y = predicted, ymin = conf.low, ymax = conf.high)) + labs (x = "出場順番" , y = "平均得点（100点満点）" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 )
各point-rangeを線で繋ぐためにはgeom_line()オブジェクトで折れ線グラフを追加する。
%>% ggplot () + geom_pointrange (aes (x = Order, y = predicted, ymin = conf.low, ymax = conf.high)) + geom_line (aes (x = Order, y = predicted)) + labs (x = "出場順番" , y = "平均得点（100点満点）" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 )
この場合、geom_pointrange()とgeom_line()はxとyに対して同じマッピングを共有するため、ggplot()内部でマッピングした方がより効率的な書き方だ。
%>% ggplot (aes (x = Order, y = predicted)) + geom_pointrange (aes (ymin = conf.low, ymax = conf.high)) + geom_line () + labs (x = "出場順番" , y = "平均得点（100点満点）" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 )
ここで予測値の図における95%信頼区間の意味について説明する。先程の図における右上がりの直線が実は「回帰直線 （regression line）」であり、2つの変数の関係を1次関数で示したのもである。しかし、1次関数の係数（\(\alpha\) や\(\beta_1\) など）には不確実性が存在し、だからこそ統計的有意性検定を行うことである。これらの係数が1次関数の切片と傾きであるが、不確実性が存在する以上、回帰直線にも不確実性が存在すると考えた方が自然だろう。この95%信頼区間の内部には以下のように無数の直線が引ける。
青色の直線が推定結果から得られた回帰直線であるが、青以外の直線も95%信頼区間内部に存在する。これらの直線は今回の推定結果から得られた直線ではないが、あり得る 回帰直線である。つまり、もう一度、2013年から同じやり方でM-1グランプリを実施し、同じやり方でデータ収集&分析をすれば、緑やオレンジ色のような回帰直線が得られたとしてもそこまで不思議ではないことを意味する。ちなみにここで最も説明変数の影響力が小さいのはオレンジ色の直線である。一つ言えるのは、説明変数の傾き係数が統計的有意であれば、いくら効果量が小さい直線であっても、右上がり関係を示すことである（つまり、右下がりにはならない）。一方、統計的に有意でない場合は、区間内に右上がりと右下がりの直線が共存することになる。
ロジスティック回帰分析
ロジスティック回帰分析は線形回帰分析とは区別される手法であるが、説明変数の値ごとの予測値を可視化する方法は同じだ。つまり、predictions()で予測値を計算し、その計算結果が格納された表形式データ（data.frame、またはtibbleと呼ばれる）と{ggplot2}を用いた作図をするだけだ。まずはOrderの値が1、2、…、10の場合のファイナルステージへの進出確率の予測値を計算し、logit_predに格納する。このlogit_predを使ってそのまま作図へ移行しても良いが、自分で具体的な数値を確認したい場合は、必要な列だけ抽出しておくと良い。
<- predictions (fit_logit, newdata = datagrid (Order = 1 : 10 ))<- logit_pred %>% select (Order, predicted, std.error, conf.low, conf.high)
Order predicted std.error conf.low conf.high
1 1 0.09377992 0.06076774 0.02484884 0.2959036
2 2 0.11999376 0.06411783 0.03982630 0.3095141
3 3 0.15230355 0.06535803 0.06245367 0.3264132
4 4 0.19142073 0.06456340 0.09462983 0.3490452
5 5 0.23776669 0.06327313 0.13596282 0.3820889
6 6 0.29129134 0.06526040 0.18113035 0.4330217
7 7 0.35131170 0.07493977 0.22135429 0.5078080
8 8 0.41643237 0.09301367 0.25206946 0.6017430
9 9 0.48460633 0.11581347 0.27479749 0.6999854
10 10 0.55335806 0.13842452 0.29245683 0.7878439
続いて、{ggplot2}を用いて作図をする。基本的に線形回帰分析と同じコードであるが、今回の縦軸は「平均得点」の予測値でなく、「ファイナルステージへの進出確率」の予測値であるため、軸ラベルを修正する。
%>% ggplot () + geom_pointrange (aes (x = Order, y = predicted, ymin = conf.low, ymax = conf.high)) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 )
交互作用
予測値（調整変数がダミー変数）
交差項を含むモデル（線形回帰・ロジスティック回帰両方）の場合でも基本的な手順はこれまでと同じである。しかし、主な説明変数（ここではOrder）が応答変数に与える影響は調整変数（ここではFirst、またはDuration）によって異なる。したがって、調整変数の値ごとに別々の回帰直線を計算する必要がある。そのためにはpredictions()のnewdata = datagrid()内に主な説明変数だけでなく、調整変数の値も指定する必要がある。たとえば、調整変数がダミー変数（ここではFirst）の場合、Firstの値が0と1の場合に分けて、Orderの値が1、2、3、…、10の場合の予測値を計算する必要がある。
<- predictions (fit_inter1, newdata = datagrid (Order = 1 : 10 ,First = 0 : 1 ))
rowid type predicted std.error statistic p.value conf.low conf.high
1 1 response 90.86470 0.8053717 112.8233 0 89.25070 92.47870
2 2 response 89.52837 0.7705319 116.1903 0 87.98419 91.07255
3 3 response 91.07488 0.6783848 134.2525 0 89.71537 92.43439
4 4 response 89.82365 0.6588679 136.3303 0 88.50325 91.14405
5 5 response 91.28506 0.5672576 160.9235 0 90.14825 92.42187
6 6 response 90.11893 0.5596164 161.0370 0 88.99743 91.24043
7 7 response 91.49524 0.4830628 189.4065 0 90.52716 92.46332
8 8 response 90.41421 0.4805313 188.1547 0 89.45121 91.37722
9 9 response 91.70542 0.4414875 207.7192 0 90.82066 92.59018
10 10 response 90.70950 0.4328119 209.5818 0 89.84212 91.57687
11 11 response 91.91560 0.4543851 202.2857 0 91.00499 92.82621
12 12 response 91.00478 0.4271044 213.0739 0 90.14884 91.86071
13 13 response 92.12578 0.5177003 177.9520 0 91.08829 93.16327
14 14 response 91.30006 0.4649583 196.3618 0 90.36826 92.23186
15 15 response 92.33596 0.6160799 149.8766 0 91.10131 93.57061
16 16 response 91.59534 0.5372434 170.4913 0 90.51868 92.67200
17 17 response 92.54614 0.7355871 125.8126 0 91.07199 94.02029
18 18 response 91.89062 0.6322585 145.3371 0 90.62355 93.15770
19 19 response 92.75632 0.8675340 106.9195 0 91.01774 94.49490
20 20 response 92.18591 0.7413146 124.3546 0 90.70028 93.67153
Score_Mean Duration Order First
1 91.32857 10.6 1 0
2 91.32857 10.6 1 1
3 91.32857 10.6 2 0
4 91.32857 10.6 2 1
5 91.32857 10.6 3 0
6 91.32857 10.6 3 1
7 91.32857 10.6 4 0
8 91.32857 10.6 4 1
9 91.32857 10.6 5 0
10 91.32857 10.6 5 1
11 91.32857 10.6 6 0
12 91.32857 10.6 6 1
13 91.32857 10.6 7 0
14 91.32857 10.6 7 1
15 91.32857 10.6 8 0
16 91.32857 10.6 8 1
17 91.32857 10.6 9 0
18 91.32857 10.6 9 1
19 91.32857 10.6 10 0
20 91.32857 10.6 10 1
作図に必要な列（+標準誤差であるstd.error列）のみを抽出する。今回の例はOrder以外にもFirstの値も2パターン（0と1）あるため、First列も抽出する。
<- inter1_pred %>% select (Order, First, predicted, std.error, conf.low, conf.high)
Order First predicted std.error conf.low conf.high
1 1 0 90.86470 0.8053717 89.25070 92.47870
2 1 1 89.52837 0.7705319 87.98419 91.07255
3 2 0 91.07488 0.6783848 89.71537 92.43439
4 2 1 89.82365 0.6588679 88.50325 91.14405
5 3 0 91.28506 0.5672576 90.14825 92.42187
6 3 1 90.11893 0.5596164 88.99743 91.24043
7 4 0 91.49524 0.4830628 90.52716 92.46332
8 4 1 90.41421 0.4805313 89.45121 91.37722
9 5 0 91.70542 0.4414875 90.82066 92.59018
10 5 1 90.70950 0.4328119 89.84212 91.57687
11 6 0 91.91560 0.4543851 91.00499 92.82621
12 6 1 91.00478 0.4271044 90.14884 91.86071
13 7 0 92.12578 0.5177003 91.08829 93.16327
14 7 1 91.30006 0.4649583 90.36826 92.23186
15 8 0 92.33596 0.6160799 91.10131 93.57061
16 8 1 91.59534 0.5372434 90.51868 92.67200
17 9 0 92.54614 0.7355871 91.07199 94.02029
18 9 1 91.89062 0.6322585 90.62355 93.15770
19 10 0 92.75632 0.8675340 91.01774 94.49490
20 10 1 92.18591 0.7413146 90.70028 93.67153
続いて、First列をリコーディングする。分析する側はFirstの値が1だと「初出場」、0だと「出場経験あり」とうことが分かるが、一般読者は分からないだろう。ここではFirstをfactor化し、0と1にそれぞれ「出場経験あり」、「初出場」とラベルを付けておこう。
<- inter1_pred %>% mutate (First = factor (First, levels = c (0 , 1 ), labels = c ("出場経験あり" , "初出場" )))
Order First predicted std.error conf.low conf.high
1 1 出場経験あり 90.86470 0.8053717 89.25070 92.47870
2 1 初出場 89.52837 0.7705319 87.98419 91.07255
3 2 出場経験あり 91.07488 0.6783848 89.71537 92.43439
4 2 初出場 89.82365 0.6588679 88.50325 91.14405
5 3 出場経験あり 91.28506 0.5672576 90.14825 92.42187
6 3 初出場 90.11893 0.5596164 88.99743 91.24043
7 4 出場経験あり 91.49524 0.4830628 90.52716 92.46332
8 4 初出場 90.41421 0.4805313 89.45121 91.37722
9 5 出場経験あり 91.70542 0.4414875 90.82066 92.59018
10 5 初出場 90.70950 0.4328119 89.84212 91.57687
11 6 出場経験あり 91.91560 0.4543851 91.00499 92.82621
12 6 初出場 91.00478 0.4271044 90.14884 91.86071
13 7 出場経験あり 92.12578 0.5177003 91.08829 93.16327
14 7 初出場 91.30006 0.4649583 90.36826 92.23186
15 8 出場経験あり 92.33596 0.6160799 91.10131 93.57061
16 8 初出場 91.59534 0.5372434 90.51868 92.67200
17 9 出場経験あり 92.54614 0.7355871 91.07199 94.02029
18 9 初出場 91.89062 0.6322585 90.62355 93.15770
19 10 出場経験あり 92.75632 0.8675340 91.01774 94.49490
20 10 初出場 92.18591 0.7413146 90.70028 93.67153
それでは作図をしてみよう。作図のコードもこれまでと大きく変わらない。しかし、今回は「初出場」か「出場経験あり」かで別のpoint-rangeが出力される。この2種類のpoint-rangeを今回は色分け してみよう。特定の変数の値に応じて色分けをする場合はaes()の内部でcolor変数にその変数を指定する。ここだと、color = Firstとなる。
%>% ggplot () + geom_pointrange (aes (x = Order, y = predicted, ymin = conf.low, ymax = conf.high,color = First)) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" , color = "" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 ) + theme (legend.position = "bottom" )
これだと予測値（点）と信頼区間（線）が重なって読みにくい。このように予測値/限界効果を2つ以上のグループに分けて可視化する場合は、線（line）と面（ribbon）を組み合わせた方が良い。まずは、折れ線グラフを追加してみよう。折れ線グラフの幾何オブジェクトはgeom_line()だ。折れ線グラフは基本的にxとyのみマッピングで良いが、ここではFirstの値ごとに線の色分けがしたいのでcolorにもマッピングする。
%>% ggplot () + geom_line (aes (x = Order, y = predicted, color = First)) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" , color = "" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 ) + theme (legend.position = "bottom" )
つづいて、信頼区間だ。信頼区間を面で示す場合はgeom_ribbon()幾何オブジェクトを使用する。マッピングはgeom_pointrange()と似ているが、yのマッピングは不要である。また、線の色はcolorにマッピングをしたが、面の色はfillである。
%>% ggplot () + geom_line (aes (x = Order, y = predicted, color = First)) + geom_ribbon (aes (x = Order, ymin = conf.low, ymax = conf.high,fill = First)) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" , color = "" , fill = "" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 ) + theme (legend.position = "bottom" )
今回は面が重なってしまう。しかも線と面が同じいろなので、線が見えなくなってしまった。この場合は面を半透明にすれば良い。すべての面に適用するものはaes()外側に書くが、透明度を調整する引数はalphaであり、0に近いほど透明になる。今回は透明度30%（alpha = 0.3）をaes()の外側に追加する。また、線の太さも太めにしたいので、geom_line()のaes()の外側にlinewidth引数を追加して線の太さを調整する。最後に、geom_line()とgeom_ribbon()の順番を変えることで、線が面より上に表示されるようにする（{ggplot2}は後で追加したレイヤーがより上に位置する）。
%>% ggplot () + geom_ribbon (aes (x = Order, ymin = conf.low, ymax = conf.high,fill = First), alpha = 0.3 ) + geom_line (aes (x = Order, y = predicted, color = First), linewidth = 1 ) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" , color = "" , fill = "" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 ) + theme (legend.position = "bottom" )
これで作図は終わったが、2つの幾何オブジェクトがxに対して同じマッピングを共有しているため、こちらはggplot()内にマッピングするとより簡潔なコードになる。
%>% ggplot (aes (x = Order)) + geom_ribbon (aes (ymin = conf.low, ymax = conf.high, fill = First), alpha = 0.3 ) + geom_line (aes (y = predicted, color = First), linewidth = 1 ) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" , color = "" , fill = "" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw (base_size = 14 ) + theme (legend.position = "bottom" )
予測値（調整変数が連続変数）
調整変数が連続変数の場合でも同じ手順で可視化する。まずは、{marginaleffects}のpredictions()で予測値を計算し、表形式オブジェクトとして格納する。ただし、Durationの値が2、3、4、…、16の場合の予測値をすべて出してしまうと15本の線（と面）が表示されるため、非常に可読性が落ちてしまう。したがって、ここではDurationの最小値（2）、最大値（16）とその中間である9に固定して予測値を計算する。計算結果はinter2_predに格納する。
<- predictions (fit_inter2, newdata = datagrid (Order = 1 : 10 ,Duration = c (2 , 9 , 16 )))
rowid type predicted std.error statistic p.value conf.low conf.high
1 1 response 88.75862 1.5824450 56.08954 0 85.58733 91.92991
2 2 response 89.86234 0.6614179 135.86319 0 88.53683 91.18785
3 3 response 90.96607 0.9336240 97.43330 0 89.09504 92.83709
4 4 response 89.18246 1.3389996 66.60380 0 86.49905 91.86588
5 5 response 90.15477 0.5613941 160.59088 0 89.02972 91.27983
6 6 response 91.12709 0.7872183 115.75834 0 89.54947 92.70471
7 7 response 89.60630 1.1184767 80.11459 0 87.36483 91.84778
8 8 response 90.44721 0.4701238 192.39020 0 89.50506 91.38936
9 9 response 91.28811 0.6636047 137.56400 0 89.95822 92.61801
10 10 response 90.03015 0.9371995 96.06295 0 88.15196 91.90834
11 11 response 90.73964 0.3937419 230.45463 0 89.95056 91.52872
12 12 response 91.44913 0.5776064 158.32432 0 90.29159 92.60668
13 13 response 90.45399 0.8215707 110.09886 0 88.80753 92.10046
14 14 response 91.03207 0.3423628 265.89363 0 90.34596 91.71818
15 15 response 91.61016 0.5472536 167.39981 0 90.51344 92.70688
16 16 response 90.87784 0.8005600 113.51783 0 89.27348 92.48219
17 17 response 91.32451 0.3279560 278.46577 0 90.66727 91.98175
18 18 response 91.77118 0.5813291 157.86443 0 90.60617 92.93619
19 19 response 91.30168 0.8809637 103.63842 0 89.53619 93.06717
20 20 response 91.61694 0.3550510 258.03880 0 90.90540 92.32848
21 21 response 91.93220 0.6700746 137.19697 0 90.58934 93.27506
22 22 response 91.72552 1.0395096 88.23923 0 89.64230 93.80875
23 23 response 91.90937 0.4156087 221.14398 0 91.07648 92.74227
24 24 response 92.09322 0.7953967 115.78276 0 90.49921 93.68724
25 25 response 92.14937 1.2467341 73.91261 0 89.65086 94.64788
26 26 response 92.20181 0.4975573 185.30891 0 91.20468 93.19893
27 27 response 92.25425 0.9428209 97.84917 0 90.36479 94.14370
28 28 response 92.57321 1.4823610 62.44984 0 89.60249 95.54393
29 29 response 92.49424 0.5920804 156.21904 0 91.30769 93.68080
30 30 response 92.41527 1.1035245 83.74556 0 90.20376 94.62678
Score_Mean First Order Duration
1 91.32857 0.5166667 1 2
2 91.32857 0.5166667 1 9
3 91.32857 0.5166667 1 16
4 91.32857 0.5166667 2 2
5 91.32857 0.5166667 2 9
6 91.32857 0.5166667 2 16
7 91.32857 0.5166667 3 2
8 91.32857 0.5166667 3 9
9 91.32857 0.5166667 3 16
10 91.32857 0.5166667 4 2
11 91.32857 0.5166667 4 9
12 91.32857 0.5166667 4 16
13 91.32857 0.5166667 5 2
14 91.32857 0.5166667 5 9
15 91.32857 0.5166667 5 16
16 91.32857 0.5166667 6 2
17 91.32857 0.5166667 6 9
18 91.32857 0.5166667 6 16
19 91.32857 0.5166667 7 2
20 91.32857 0.5166667 7 9
21 91.32857 0.5166667 7 16
22 91.32857 0.5166667 8 2
23 91.32857 0.5166667 8 9
24 91.32857 0.5166667 8 16
25 91.32857 0.5166667 9 2
26 91.32857 0.5166667 9 9
27 91.32857 0.5166667 9 16
28 91.32857 0.5166667 10 2
29 91.32857 0.5166667 10 9
30 91.32857 0.5166667 10 16
<- inter2_pred %>% select (Order, Duration, predicted, std.error, conf.low, conf.high)
Order Duration predicted std.error conf.low conf.high
1 1 2 88.75862 1.5824450 85.58733 91.92991
2 1 9 89.86234 0.6614179 88.53683 91.18785
3 1 16 90.96607 0.9336240 89.09504 92.83709
4 2 2 89.18246 1.3389996 86.49905 91.86588
5 2 9 90.15477 0.5613941 89.02972 91.27983
6 2 16 91.12709 0.7872183 89.54947 92.70471
7 3 2 89.60630 1.1184767 87.36483 91.84778
8 3 9 90.44721 0.4701238 89.50506 91.38936
9 3 16 91.28811 0.6636047 89.95822 92.61801
10 4 2 90.03015 0.9371995 88.15196 91.90834
11 4 9 90.73964 0.3937419 89.95056 91.52872
12 4 16 91.44913 0.5776064 90.29159 92.60668
13 5 2 90.45399 0.8215707 88.80753 92.10046
14 5 9 91.03207 0.3423628 90.34596 91.71818
15 5 16 91.61016 0.5472536 90.51344 92.70688
16 6 2 90.87784 0.8005600 89.27348 92.48219
17 6 9 91.32451 0.3279560 90.66727 91.98175
18 6 16 91.77118 0.5813291 90.60617 92.93619
19 7 2 91.30168 0.8809637 89.53619 93.06717
20 7 9 91.61694 0.3550510 90.90540 92.32848
21 7 16 91.93220 0.6700746 90.58934 93.27506
22 8 2 91.72552 1.0395096 89.64230 93.80875
23 8 9 91.90937 0.4156087 91.07648 92.74227
24 8 16 92.09322 0.7953967 90.49921 93.68724
25 9 2 92.14937 1.2467341 89.65086 94.64788
26 9 9 92.20181 0.4975573 91.20468 93.19893
27 9 16 92.25425 0.9428209 90.36479 94.14370
28 10 2 92.57321 1.4823610 89.60249 95.54393
29 10 9 92.49424 0.5920804 91.30769 93.68080
30 10 16 92.41527 1.1035245 90.20376 94.62678
続いて、調整変数をfactor化する。inter2_predのDuration列は2、9、16の値で構成されており（これはlevelsで指定する）、それぞれ「2年目」、「9年目」、「16年目」とラベルを付ける（これはlabelsで指定する）。
<- inter2_pred %>% mutate (Duration = factor (Duration, levels = c (2 , 9 , 16 ), labels = c ("2年目" , "9年目" , "16年目" )))
Order Duration predicted std.error conf.low conf.high
1 1 2年目 88.75862 1.5824450 85.58733 91.92991
2 1 9年目 89.86234 0.6614179 88.53683 91.18785
3 1 16年目 90.96607 0.9336240 89.09504 92.83709
4 2 2年目 89.18246 1.3389996 86.49905 91.86588
5 2 9年目 90.15477 0.5613941 89.02972 91.27983
6 2 16年目 91.12709 0.7872183 89.54947 92.70471
7 3 2年目 89.60630 1.1184767 87.36483 91.84778
8 3 9年目 90.44721 0.4701238 89.50506 91.38936
9 3 16年目 91.28811 0.6636047 89.95822 92.61801
10 4 2年目 90.03015 0.9371995 88.15196 91.90834
11 4 9年目 90.73964 0.3937419 89.95056 91.52872
12 4 16年目 91.44913 0.5776064 90.29159 92.60668
13 5 2年目 90.45399 0.8215707 88.80753 92.10046
14 5 9年目 91.03207 0.3423628 90.34596 91.71818
15 5 16年目 91.61016 0.5472536 90.51344 92.70688
16 6 2年目 90.87784 0.8005600 89.27348 92.48219
17 6 9年目 91.32451 0.3279560 90.66727 91.98175
18 6 16年目 91.77118 0.5813291 90.60617 92.93619
19 7 2年目 91.30168 0.8809637 89.53619 93.06717
20 7 9年目 91.61694 0.3550510 90.90540 92.32848
21 7 16年目 91.93220 0.6700746 90.58934 93.27506
22 8 2年目 91.72552 1.0395096 89.64230 93.80875
23 8 9年目 91.90937 0.4156087 91.07648 92.74227
24 8 16年目 92.09322 0.7953967 90.49921 93.68724
25 9 2年目 92.14937 1.2467341 89.65086 94.64788
26 9 9年目 92.20181 0.4975573 91.20468 93.19893
27 9 16年目 92.25425 0.9428209 90.36479 94.14370
28 10 2年目 92.57321 1.4823610 89.60249 95.54393
29 10 9年目 92.49424 0.5920804 91.30769 93.68080
30 10 16年目 92.41527 1.1035245 90.20376 94.62678
最後に{ggplot2}で作図する。使用する幾何オブジェクトはこれまでと同様、geom_line()とgeom_ribbon()だ。
%>% ggplot (aes (x = Order)) + geom_ribbon (aes (ymin = conf.low, ymax = conf.high, fill = Duration), alpha = 0.3 ) + geom_line (aes (y = predicted, color = Duration), linewidth = 1 ) + labs (x = "出場順番" , y = "ファイナルステージへの進出確率" , color = "結成からの経過年数" , fill = "結成からの経過年数" ) + scale_x_continuous (breaks = 1 : 10 , labels = 1 : 10 ) + theme_bw () + theme (legend.position = "bottom" )
限界効果（調整変数がダミー変数）
以上の例は{marginaleffects}パッケージのpredictions()関数を使って予測値 を計算してきたが、ここからは同じパッケージのmarginaleffects()関数を使って限界効果 を計算し、{ggplot2}で可視化する。marginaleffects()の場合、限界効果を計算した変数はvariables引数に、調整変数はこれまでと同様、newdata = datagrid()内に指定 する。限界効果を計算したい変数、つまり交差項の中で調整変数ではない ものが主な制つ名変数であり、今回の例だと出場順番（Order）である。この変数を"を囲み、variablesに割り当てる。調整変数は初出場ダミー（First）が調整変数であれば、newdata = datagrid()のdatagrid()内にFirst = 0:1を割り当てる。計算結果はinter1_ameに格納する。
<- marginaleffects (fit_inter1, variables = "Order" ,newdata = datagrid (First = 0 : 1 ))
rowid type term dydx std.error statistic p.value conf.low
1 1 response Order 0.2101801 0.1581030 1.329387 0.1837203 -0.09969606
2 2 response Order 0.2952823 0.1390614 2.123396 0.0337207 0.02272706
conf.high predicted predicted_hi predicted_lo Score_Mean Order Duration First
1 0.5200562 91.81051 91.8107 91.81051 91.32857 5.5 10.6 0
2 0.5678376 90.85714 90.8574 90.85714 91.32857 5.5 10.6 1
eps
1 9e-04
2 9e-04
作図に必要な列（と標準誤差を意味するstd.error列）を抽出する。一つ注意すべき点は予測値はpredicted列であったのに対し、限界効果の計算結果には当該列が存在せず、dydx列に格納されている点だ。また、限界効果を可視化する際、Orderの値は不要であるため、省略しても良い。
<- inter1_ame %>% select (First, dydx, std.error, conf.low, conf.high)
First dydx std.error conf.low conf.high
1 0 0.2101801 0.1581030 -0.09969606 0.5200562
2 1 0.2952823 0.1390614 0.02272706 0.5678376
続いて、inter1_ameのFirst列をfactor化する。
<- inter1_ame %>% mutate (First = factor (First, levels = 0 : 1 , labels = c ("出場経験あり" , "初出場" )))
First dydx std.error conf.low conf.high
1 出場経験あり 0.2101801 0.1581030 -0.09969606 0.5200562
2 初出場 0.2952823 0.1390614 0.02272706 0.5678376
それでは作図してみよう。調整変数が連続変数であれ、ダミー変数であれ、取り得る値が数個程度ならgeom_pointrange()を使用する。横軸（x）は調整変数、縦軸（y）は限界効果にし、限界効果の95%信頼区間の下限と上限はそれぞれyminとymaxにマッピングする。また、限界効果の統計的有意性を確認するために、y = 0の箇所に水平線を引く。水平線の幾何オブジェクトはgeom_hline()であり、()内にyintercept = 0を入れると、y = 0のところに水平線が表示される。
%>% ggplot () + geom_hline (yintercept = 0 ) + geom_pointrange (aes (x = First, y = dydx, ymin = conf.low, ymax = conf.high)) + labs (x = "出場経験" , y = "出場順番の平均限界効果" ) + theme_bw (base_size = 12 )
95％信頼区間内に0が含まれる場合、当該限界効果は統計的に有意でないことを意味する。たとえば、上の図だと初出場の場合のみ、出場順番の限界効果が統計的有意であることが分かる。つまり、「初出場の人の場合、出場順番が遅いほど平均得点が高い傾向がある」と解釈できる。一方、既に出場経験がある場合、出場順番の限界効果が統計的有意ではない。これは、「M-1グランプリの出場経験がある場合、出場順番が遅いほど平均得点が高い傾向があるとは言えない 」と解釈できる。「〜する傾向はない」ではなく、「〜する傾向があるとは言えない」と解釈する点に注意すること。一般的な統計的有意性検定において「ない」と言い切ることは至難の業である。
限界効果（調整変数が連続変数）
最後に調整変数が連続変数（ここではDuration）の場合について考えてみよう。基本的なやり方はこれまでと同じである。しかし、予測値の可視化の際、Durationの値が2、9、16の場合の予測値を計算した。なぜなら、この値が3つでなく、5つになると線の数が増えてしまい、図の可読性が低下するからだ。情報量が多少減っても、可読性の低い図は存在意義がないため、線の本数を減らしたわけだ。しかし、今回の例だと、Durationが線の色分けでなく、横軸として使用される。したがって、Durationのすべての値を対象に限界効果を計算した方が良い。datagrid()内にはDuration = c(1, 2, 3, 4, 5)と書いても良いが、Duration = 2:16の方がより効率的な書き方である。
<- marginaleffects (fit_inter2, variables = "Order" ,newdata = datagrid (Duration = 2 : 16 ))
rowid type term dydx std.error statistic p.value conf.low
1 1 response Order 0.4238440 0.2909318 1.4568503 0.14515771 -0.146371757
2 2 response Order 0.4050711 0.2624520 1.5434099 0.12273135 -0.109325467
3 3 response Order 0.3862981 0.2345512 1.6469670 0.09956483 -0.073413823
4 4 response Order 0.3675251 0.2074630 1.7715214 0.07647404 -0.039094864
5 5 response Order 0.3487522 0.1815515 1.9209546 0.05473743 -0.007082209
6 6 response Order 0.3299792 0.1573989 2.0964517 0.03604214 0.021483018
7 7 response Order 0.3112063 0.1359460 2.2891910 0.02206826 0.044757079
8 8 response Order 0.2924333 0.1186661 2.4643379 0.01372666 0.059852078
9 9 response Order 0.2736604 0.1075889 2.5435731 0.01097251 0.062789886
10 10 response Order 0.2548874 0.1047022 2.4344034 0.01491636 0.049674854
11 11 response Order 0.2361144 0.1106487 2.1339099 0.03285016 0.019246897
12 12 response Order 0.2173415 0.1241659 1.7504126 0.08004714 -0.026019124
13 13 response Order 0.1985685 0.1431244 1.3873840 0.16532471 -0.081950181
14 14 response Order 0.1797956 0.1656667 1.0852847 0.27779558 -0.144905251
15 15 response Order 0.1610226 0.1905249 0.8451524 0.39802578 -0.212399410
conf.high predicted predicted_hi predicted_lo Score_Mean Order First
1 0.9940598 90.66591 90.66630 90.66591 91.32857 5.5 0.5166667
2 0.9194676 90.73911 90.73948 90.73911 91.32857 5.5 0.5166667
3 0.8460100 90.81231 90.81266 90.81231 91.32857 5.5 0.5166667
4 0.7741452 90.88550 90.88584 90.88550 91.32857 5.5 0.5166667
5 0.7045866 90.95870 90.95901 90.95870 91.32857 5.5 0.5166667
6 0.6384754 91.03190 91.03219 91.03190 91.32857 5.5 0.5166667
7 0.5776555 91.10509 91.10537 91.10509 91.32857 5.5 0.5166667
8 0.5250145 91.17829 91.17855 91.17829 91.32857 5.5 0.5166667
9 0.4845308 91.25149 91.25173 91.25149 91.32857 5.5 0.5166667
10 0.4600999 91.32468 91.32491 91.32468 91.32857 5.5 0.5166667
11 0.4529820 91.39788 91.39809 91.39788 91.32857 5.5 0.5166667
12 0.4607021 91.47108 91.47127 91.47108 91.32857 5.5 0.5166667
13 0.4790872 91.54427 91.54445 91.54427 91.32857 5.5 0.5166667
14 0.5044964 91.61747 91.61763 91.61747 91.32857 5.5 0.5166667
15 0.5344446 91.69067 91.69081 91.69067 91.32857 5.5 0.5166667
Duration eps
1 2 9e-04
2 3 9e-04
3 4 9e-04
4 5 9e-04
5 6 9e-04
6 7 9e-04
7 8 9e-04
8 9 9e-04
9 10 9e-04
10 11 9e-04
11 12 9e-04
12 13 9e-04
13 14 9e-04
14 15 9e-04
15 16 9e-04
続いて、作図に必要な列を出力する。今回の場合、連続変数であるため、factor化は不要 だ。
<- inter2_ame %>% select (Duration, dydx, std.error, conf.low, conf.high)
Duration dydx std.error conf.low conf.high
1 2 0.4238440 0.2909318 -0.146371757 0.9940598
2 3 0.4050711 0.2624520 -0.109325467 0.9194676
3 4 0.3862981 0.2345512 -0.073413823 0.8460100
4 5 0.3675251 0.2074630 -0.039094864 0.7741452
5 6 0.3487522 0.1815515 -0.007082209 0.7045866
6 7 0.3299792 0.1573989 0.021483018 0.6384754
7 8 0.3112063 0.1359460 0.044757079 0.5776555
8 9 0.2924333 0.1186661 0.059852078 0.5250145
9 10 0.2736604 0.1075889 0.062789886 0.4845308
10 11 0.2548874 0.1047022 0.049674854 0.4600999
11 12 0.2361144 0.1106487 0.019246897 0.4529820
12 13 0.2173415 0.1241659 -0.026019124 0.4607021
13 14 0.1985685 0.1431244 -0.081950181 0.4790872
14 15 0.1797956 0.1656667 -0.144905251 0.5044964
15 16 0.1610226 0.1905249 -0.212399410 0.5344446
続いて作図をする。今回は調整変数が連続変数であるため、線（geom_line()）と面（geom_ribbon()）の図を出力してみよう。面の透明度指定を忘れないこと（一度、alpha = 0.3なしでコードを走らせてみても良いだろう）。
%>% ggplot (aes (x = Duration)) + geom_hline (yintercept = 0 ) + geom_ribbon (aes (y = dydx, ymin = conf.low, ymax = conf.high),alpha = 0.3 ) + geom_line (aes (y = dydx), linewidth = 1 ) + labs (x = "結成からの経過年数" , y = "出場順番の平均限界効果" ) + theme_bw (base_size = 14 )
今回、Durationは2から16まで、つまり14の値で構成されている。この場合、pointrangeプロットを作成しても良いかも知れない。むしろ、結果を解釈する時にはpointrangeプロットの方がやりやすい。
%>% ggplot () + geom_hline (yintercept = 0 ) + geom_pointrange (aes (x = Duration, y = dydx, ymin = conf.low, ymax = conf.high)) + labs (x = "結成からの経過年数" , y = "出場順番の平均限界効果" ) + theme_bw (base_size = 14 )
横軸の目盛りを1年刻みに修正してみよう。
%>% ggplot () + geom_hline (yintercept = 0 ) + geom_pointrange (aes (x = Duration, y = dydx, ymin = conf.low, ymax = conf.high)) + labs (x = "結成からの経過年数" , y = "出場順番の平均限界効果" ) + scale_x_continuous (breaks = 2 : 16 , labels = 2 : 16 ) + theme_bw (base_size = 14 )
限界効果が統計的に有意な箇所はDurationの値が7以上、12以下の時である。つまり、「コンビー結成から7年以上、12年以下経過した場合、出場順番が遅いほど平均得点が高い傾向がある」と解釈できる。また、「コンビー結成から6年以下、または13年以上経過した場合、出場順番が遅いほど平均得点が高い傾向があるとは言えない」といった解釈もできる。
おまけに統計的有意性検定の結果も同時に示す方法も紹介する。具体的には限界効果が統計的有意である場合は黒に、非有意である場合はグレーに出力する。コードの説明は割愛するため、興味のある方はぜひコードの中身を解読してみよう。
%>% mutate (Sig = if_else (conf.low * conf.high > 0 , "Sig" , "Insig" )) %>% ggplot () + geom_hline (yintercept = 0 ) + geom_pointrange (aes (x = Duration, y = dydx, ymin = conf.low, ymax = conf.high,color = Sig)) + labs (x = "結成からの経過年数" , y = "出場順番の平均限界効果" ) + scale_x_continuous (breaks = 2 : 16 , labels = 2 : 16 ) + scale_color_manual (values = c ("Sig" = "black" , "Insig" = "gray70" )) + theme_bw (base_size = 14 ) + theme (legend.position = "none" )
