第13回講義資料

回帰分析の可視化

スライド

新しいタブで開く

セットアップ

 本日使用するパッケージとLMSからダウンロードして実習用データを読み込む。読み込むデータはdfに格納し、中身を確認してみよう。

Code 01
library(tidyverse)
library(summarytools)
library(marginaleffects)

df <- read_csv("Data/M1.csv")

df
# 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)列の記述統計は不要であるため、予め除外しておこう。

Code 02
df %>%
  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: 線形回帰分析 + 交互作用(連続変数)
    • チョ制変数は「コンビー結成からの経過年数」
Code 03
fit_lm     <- lm(Score_Mean ~ Order + First + Duration, data = df)
fit_logit  <- glm(Final ~ Order + First + Duration, data = df, 
                  family = binomial("logit"))
fit_inter1 <- lm(Score_Mean ~ Order * First + Duration, data = df)
fit_inter2 <- lm(Score_Mean ~ Order * Duration + First, data = df)

 {modelsummary}のmodelsummary()関数を使って、4つのモデルを概観してみよう。

Code 04
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引数はnewdatadatagrid()関数を割り当てたものである。datagrid()の中には説明変数名 = 値を指定する。1のみならOrder = 1の場合の予測値が計算され、c(1, 2, 3)ならOrderの値が1、2、3の場合の予測値が計算される。公差1の等差数列は:演算子で表現することもできるので、今回は1:10とする。計算結果はlm_predに格納し、出力する。

Code 05
lm_pred <- predictions(fit_lm, newdata = datagrid(Order = 1:10))

lm_pred
   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を上書きしよう。

Code 06
lm_pred <- lm_pred %>%
  select(Order, predicted, std.error, conf.low, conf.high)

lm_pred
   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()内部でマッピングを行う。最低でもxyyminymaxのマッピングが必要であり、xyはそれぞれ横軸と縦軸に対応する変数、yminymaxは区間の下限と上限に対応する変数を指定する。また、横軸の目盛りを1刻みにするため、scale_x_continuous()レイターを足す。scale_x_continuous()breaksに目盛りの位置を、labelsには各目盛りに付けるラベルのベクトルを割り当てる。breakslabelsは同じ長さのベクトルである必要がある。

Code 07
lm_pred %>%
  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()オブジェクトで折れ線グラフを追加する。

Code 08
lm_pred %>%
  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()xyに対して同じマッピングを共有するため、ggplot()内部でマッピングした方がより効率的な書き方だ。

Code 09
lm_pred %>%
  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を使ってそのまま作図へ移行しても良いが、自分で具体的な数値を確認したい場合は、必要な列だけ抽出しておくと良い。

Code 10
logit_pred <- predictions(fit_logit, newdata = datagrid(Order = 1:10))

logit_pred <- logit_pred %>%
  select(Order, predicted, std.error, conf.low, conf.high)

logit_pred
   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}を用いて作図をする。基本的に線形回帰分析と同じコードであるが、今回の縦軸は「平均得点」の予測値でなく、「ファイナルステージへの進出確率」の予測値であるため、軸ラベルを修正する。

Code 10
logit_pred %>%
  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の場合の予測値を計算する必要がある。

Code 11
inter1_pred <- predictions(fit_inter1, newdata = datagrid(Order = 1:10,
                                                          First = 0:1))

inter1_pred
   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列も抽出する。

Code 12
inter1_pred <- inter1_pred %>%
  select(Order, First, predicted, std.error, conf.low, conf.high)

inter1_pred
   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にそれぞれ「出場経験あり」、「初出場」とラベルを付けておこう。

Code 13
inter1_pred <- inter1_pred %>%
  mutate(First = factor(First, levels = c(0, 1), labels = c("出場経験あり", "初出場")))

inter1_pred
   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となる。

Code 14
inter1_pred %>%
  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()だ。折れ線グラフは基本的にxyのみマッピングで良いが、ここではFirstの値ごとに線の色分けがしたいのでcolorにもマッピングする。

Code 15
inter1_pred %>%
  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である。

Code 16
inter1_pred %>%
  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}は後で追加したレイヤーがより上に位置する)。

Code 17
inter1_pred %>%
  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()内にマッピングするとより簡潔なコードになる。

Code 18
inter1_pred %>%
  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に格納する。

Code 19
inter2_pred <- predictions(fit_inter2, newdata = datagrid(Order = 1:10,
                                                          Duration = c(2, 9, 16)))

inter2_pred
   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
Code 20
inter2_pred <- inter2_pred %>%
  select(Order, Duration, predicted, std.error, conf.low, conf.high)

inter2_pred
   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_predDuration列は2、9、16の値で構成されており(これはlevelsで指定する)、それぞれ「2年目」、「9年目」、「16年目」とラベルを付ける(これはlabelsで指定する)。

Code 21
inter2_pred <- inter2_pred %>%
  mutate(Duration = factor(Duration, levels = c(2, 9, 16), 
                           labels = c("2年目", "9年目", "16年目")))

inter2_pred
   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()だ。

Code 22
inter2_pred %>%
  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に格納する。

Code 23
inter1_ame <- marginaleffects(fit_inter1, variables = "Order",
                              newdata = datagrid(First = 0:1))

inter1_ame
  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の値は不要であるため、省略しても良い。

Code 24
inter1_ame <- inter1_ame %>%
  select(First, dydx, std.error, conf.low, conf.high)

inter1_ame
  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_ameFirst列をfactor化する。

Code 25
inter1_ame <- inter1_ame %>%
  mutate(First = factor(First, levels = 0:1, labels = c("出場経験あり", "初出場")))

inter1_ame
         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%信頼区間の下限と上限はそれぞれyminymaxにマッピングする。また、限界効果の統計的有意性を確認するために、y = 0の箇所に水平線を引く。水平線の幾何オブジェクトはgeom_hline()であり、()内にyintercept = 0を入れると、y = 0のところに水平線が表示される。

Code 26
inter1_ame %>%
  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の方がより効率的な書き方である。

Code 27
inter2_ame <- marginaleffects(fit_inter2, variables = "Order",
                              newdata = datagrid(Duration = 2:16))

inter2_ame
   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化は不要だ。

Code 28
inter2_ame <- inter2_ame %>%
  select(Duration, dydx, std.error, conf.low, conf.high)

inter2_ame
   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なしでコードを走らせてみても良いだろう)。

Code 29
inter2_ame %>%
  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プロットの方がやりやすい。

Code 30
inter2_ame %>%
  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年刻みに修正してみよう。

Code 30
inter2_ame %>%
  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年以上経過した場合、出場順番が遅いほど平均得点が高い傾向があるとは言えない」といった解釈もできる。

 おまけに統計的有意性検定の結果も同時に示す方法も紹介する。具体的には限界効果が統計的有意である場合は黒に、非有意である場合はグレーに出力する。コードの説明は割愛するため、興味のある方はぜひコードの中身を解読してみよう。

Code 31
inter2_ame %>%
  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")