StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

ggplot2で2軸グラフを描く時の軸スケーリングの作業

ggplot2で2軸のグラフを描くときは、先日のエントリでも書いたように、ggplot2自身は左軸(第1軸)と右軸(第2軸)を別々の情報として持つことはできないので、左軸と右軸の尺度の違いを自分で設定して変換しなければならない。
あとで使いまわすので、このスケーリングの作業を行うスケーラを以下のように定義しておいた。
最初は1行で書いてggplotの描画の中に埋め込んでたけどあとで混乱しないように分けて書いておいた。

library(ggplot2)
data(airquality)  # 練習用データのよみこみ

### 変数を追加
# x軸用に月と日の列を日付の変数にしておく(2行に分けてかいてる)
airquality <- airquality %>%
  mutate(Date = paste(as.character(Month), '/', as.character(Day), sep='')) %>%
  mutate(Date = as.Date(Date, format='%m/%d'))

### y1, y2の目盛り範囲を決めておく
# ここでは恣意的に決めてるが、最大と最小を取るとかでもいいと思う
y1.lim <- c(0, 25)
y2.lim <- c(50, 100)

### スケーラの関数を書いておく
# 上でつくった、y1とy2の目盛りの範囲を定めた要素数2のベクトルを使って、
# いい感じにスケールを合わせる。

# 変数のスケーラ。
# pにy2の値ベクトルを与えると、y1の尺に合わせた数字に変換。
# y2をまずゼロ基準に戻し、y2とy1のlimの幅の比でスケーリング
# した後で、y1のゼロ基準からの乖離分を足す。

variable_scaler <- function(p, lim1, lim2){
  to_zero <- p-lim2[1]
  y1_range <- lim1[2]-lim1[1]
  y2_range <- lim2[2]-lim2[1]
  scaled <- to_zero*y1_range/y2_range
  from_zero <- scaled + lim1[1]
  return(from_zero)
}

# 第2軸の目盛りのスケーラ。
# pは、sec_axis()の'.'になる。y1の目盛りをy2の目盛りに読み替えるもの。
# y1の目盛りをまずゼロ基準に戻し、y1とy2のlimの幅の比スケーリング
# した後で、y2の目盛りのゼロ基準からの乖離分を足す。

axis_scaler <- function(p, lim1, lim2){
  to_zero <- p-lim1[1]
  y1_range <- lim1[2]-lim1[1]
  y2_range <- lim2[2]-lim2[1]
  scaled <- to_zero*y2_range/y1_range
  from_zero <- scaled + lim2[1]
  return(from_zero)
}

### 描画してみる
airquality %>%
  ggplot(aes(x=Date)) +
  geom_line(aes(y=Wind, colour='Wind')) + 
  geom_line(aes(y=variable_scaler(Temp, y1.lim, y2.lim), colour='Temp')) + 
  scale_y_continuous(limit=y1.lim,  # 第1軸の範囲
                     breaks=c(0, 5, 10, 15, 20, 25),  # 第1軸の目盛り
                     sec.axis=sec_axis(
                       ~(axis_scaler(., y1.lim, y2.lim)), # 軸スケーリング
                       breaks=c(50, 60, 70, 80, 90,100), # 第2軸の目盛り
                       name='Temp')  # y2のラベルはここで設定する
                     ) + 
  labs(title = 'DUAL AXIS CHART', 
       x='Date',
       y='Wind',
       colour = 'Variable')


もちろんy2のラベルもscale_y_continuousの中でnameをつかって指定してもよい。
2つ目の折れ線のyを指定するときに変数スケーラを使い、sec.axisを指定するときにformulaの中で軸スケーラを使う。


f:id:midnightseminar:20200908181313p:plain

折れ線グラフの端っこにラベルを付けるやつ(ggplot2)

f:id:midnightseminar:20200831132131p:plain


最近になって遅ればせながらggplot2を頻繁に使うようになってきました。
で、↑こういうふうに、折れ線の端っこにラベルを置きたいと思いました。白黒の記事原稿で4本もの折れ線を重ねるのは見づらいのでそもそもやめたほうがいいですが、人生いろいろあるわけです。
そして、ggplot2の線種は見分けがつきにくいので、凡例だけで示すのは難しい。


ググると解説ブログがいくつか見つかりますが、パット見では何をやっているのか意味がわかりやすくない気もしたので、メモしておきます。
考え方としては、

  • ggplot2の他にggrepelパッケージを入れておく。これは互いに重ならない「いい感じのラベル」を書くときに使うもの。
  • x軸の範囲を左右に少し広げておく。
  • 折れ線の末端を「点」とみなして、そこに新たに散布図を描くようなイメージで、geom_text_repelでテキストを置く。
  • 散布図にnudge(点とラベルの距離)を設定することで、強制的にラベルを左右方向にズラす。するとggrepelの機能で、「点」とラベルが線で結ばれるようになる。(点とラベルというか、正確に言うと、geom_text_repelは座標に文字を配置する散布図を描くものだが、座標で指定した位置が点として認識されている。)


という手順になります。「折れ線の端を点とみなした散布図を描いて、そこに指示線つきのラベルを添える」という考え方がポイントだと思います。
イメージとしては、


f:id:midnightseminar:20200831133007p:plain


こういうタイプのグラフをまず思い浮かべればよい。この、ラベルと点が線で繋がれてる感じのものを、折れ線グラフの端っこに重ねてやって、派手な赤点を目立たないようにしてやればいいわけです。


冒頭のグラフは以下のようなコードで描きましたが、

library(ggplot2)
library(ggrepel)

dat %>%
  ggplot(aes(x=YEAR, y=FREQ)) +
  geom_line(aes(size=CODE, linetype=CODE, colour=CODE)) +   # ここ書いとかないとmanual設定も動かない
  theme_classic() +
  geom_text_repel(
    data = dat %>% filter(YEAR==max(YEAR)),    # 折れ線の右端にうつ点
    aes(x=YEAR, y=FREQ, label = CODE),
    nudge_x = 4,
    segment.alpha = 0.5,
    size = 6,
    family="MS Gothic") + 
  geom_text_repel(
    data = dat %>% filter(YEAR==min(YEAR)),
    aes(x=YEAR, y=FREQ, label = CODE),
    nudge_x = -4,                 #  点とラベルの距離
    segment.alpha = 0.5,     # 間に引かれる線をすこし薄くする
    size = 6,                          # ラベルの大きさ
    family="MS Gothic") +  # 日本語を表示するときはフォント指定しておく
  scale_linetype_manual(values = c("solid", "dashed", "solid", "twodash")) +         # 線のタイプ
  scale_size_manual(values = c(1.1,0.9,1.3,0.9)) +                                                      # 線の太さ
  scale_color_manual(values = c('#000000','#000000','#BBBBBB','#000000')) +  # 線の色
  scale_x_continuous(limits=c(1960, 2023),
                     breaks=c(1970, 1980, 1990, 2000, 2010, 2020)) + 
  theme(text = element_text(family="MS Gothic"),
        axis.text=element_text(size=15,color="black"),
        axis.title.y = element_text(size=15,color="black"),
        legend.position = 'none') + 
  labs(x="", y="\n頻度/掲載論文数\n", color = "")
  • geom_lineの行で、線種、色、太さを分けるグループを指定しておく。
  • theme_classicはシンプルなスタイルを選択してる設定。
  • geom_text_repelを2回やってますが、1つ目はxが最大の点を利用して右側にラベルを書く作業、2つ目はxが最小の点を利用して左側にラベルを書く作業です。要するにこのグラフは、1つの折れ線と2つの散布図が重なったものってことです。
  • nudgeしない場合、点の座標にそのままラベルが書かれることになり、ラベルのテキスト同士が重なる場合だけggrepelの機能でずらされて、場合によっては指示線が出るのですが、正負方向に4年分だけnudgeする設定にしておくことで、強制的にこの指示線を出すようにします。
  • scale_linetype_manual、scale_size_manual、scale_color_manualは、線の種類・太さ・色を手動設定するもので、これを設定しなければ、グループ変数(ここではCODE)ごとに適当に割当られます。
  • x軸のデータは1965年から2017年までしかないのですが、ラベルを表示する領域を確保するために、scale_x_continuousのlimitsで左右の範囲を少し広げています。
  • legend.position = 'none'で凡例は無しにしている。


f:id:midnightseminar:20200831132131p:plain

Rで棒グラフと折れ線グラフを重ねた2軸グラフを描く

さっき、Rで棒グラフと折れ線グラフを重ねたものを作ろうとして、けっこう手間取りました。最終的に描いたのは↓のようなものなのですが。


f:id:midnightseminar:20200811212316p:plain


「2軸グラフの書き方」「種類の異なるグラフの重ね方」についていろいろ調べたところ、barplot()とplot()を組み合わせるやり方もあるんですが、ggplot2でやるほうがやりやすかったです。以下、まず単なる2軸グラフの書き方をおさらいした後で、棒と折れ線を組み合わせる方法をメモしておきます。
 
 

2軸グラフの作り方

左右の軸をつかってスケールの異なるグラフを重ねたいだけなら、さほど難しくはなく、plot()ですぐできます。
例えば折れ線グラフを2つ重ねたいのだとしたら、まず1つめの変数のグラフを軸無しで書いたあとに、x軸と左軸を描く。その後、par(new=T)で、2つめの変数のグラフを軸なしで重ね描きした上で、右軸を描けばいいです。
で、最後にbox()で枠線を入れ、凡例を付けたければ付けます。

以下、適当に乱数でつくった変数で実行例を書いておきますが、冒頭の成果物にあわせてx軸の変数を日付型にしてるので、そこだけ多少ややこしくなっています。ただの数字であれば、x軸を描くときにaxis.Dateを使う必要はないです。

# 練習用の変数を適当に乱数でこしらえる
x.date <- seq(from=as.Date('2020-04-01'), to=as.Date('2020-05-31'), by=1)
y1 <- sin(seq(length(x))/7)+rnorm(n=length(x), mean=0, sd=0.2) + 2
y2 <- sin(seq(length(x))/5)*100 + rnorm(n=length(x), mean=0, sd=20) + 200

# y軸の範囲をそれぞれ決めておく
y1.lim <- c(min(y1), max(y1))
y2.lim <- c(min(y2), max(y2))

# グラフの左右の余白を少し多めにするためparを設定しとく(特に右)
par(oma = c(0, 1, 0, 3))

# 1枚目のグラフをかく(y1)
# いったん軸なしにするためaxes=Fにしてる
plot(x=x.date, y=y1, ylim=y1.lim, type='l', lwd=1.5,
     xlab='Date', ylab='y1', 
     axes = F,
     main='plot関数での二軸グラフ')

# x軸を追加(日付データなのでaxis.Dateを使う)
axis.Date(1,at=seq(min(x.date), max(x.date),"week"),format="%m/%d")

# 左の軸をかく
axis(2)

# 2枚目のグラフを重ねる(y2)
par(new=T)
plot(x=x, y=y2, ylim=y2.lim, type='l', lty='dotted', lwd=1.5,
     xlab='', ylab='', axes = F)

# 右側の軸の名前をかく
mtext('y2', side = 4, line = 3)

# 右側の軸を表示
axis(4)

# 枠をかく
box()

# 凡例
legend("bottomleft", legend = c("y1", "y2"), lty = c('solid', 'dotted'), lwd=1.5)


f:id:midnightseminar:20200811212338p:plain
 
 

ggplot2で棒グラフを折れ線グラフを重ねる

さて、今度は種類の異なるグラフを重ねるやり方ですが、barplot()とplot()を重ねるやり方だと、今回は日付のデータを使ってることもあって、x軸のコントロールが難しかったので、ggplot2でやることにしました。
あと、今日気づきましたが、MacのRStudioだとplotするときに余白が足りませんというエラーが出まくるのが、ggplot2だと出ないんですね。いままで基本的にplot派でしたが、ggplot派に改宗しようかなと思いました……。


さて作図ですが、半分ぐらいは、西浦博(8割おじさん)氏が5月ぐらいに出していた、新型コロナの実行再生算数を計算するプログラムの作図のところを参考にさせて頂きました。
最大のポイントは、左右の軸のスケールの調整です。
上述の「plotの2軸化」の場合、右側の軸の幅は、2つめの変数の値の幅がそのまま反映されていました。plotはそもそもスケールの異なるグラフを重ねることができるようになってて*1、今回の場合でいうと、後で描いたグラフの主軸を単に右側表示にしただけというわけです。


一方、ggplot2の場合は、2軸グラフを描くときも、y軸のスケール(軸の最小値と最大値)はあくまで共通になります。今回の場合、y軸の縦幅はあくまで、1つめの変数y1の最小〜最大の幅に合わせた尺度で固定される感じになります。
じゃあ、そこにどうやって2つ目の変数を重ねるのかというと、両変数の縮尺を先に計算しておいて、2つ目の変数を1つ目の変数にあわせて縮めたり伸ばしたりして収めるわけです。
で、その後で、右側の軸のところに好きなように(ただし左軸からの変換という形で)目盛りを打つことができるので、この目盛りを、もとの第2変数に対応するものにしておけばよい。


少し分かりづらいですが、たとえば左の軸で表現したい第1変数が0〜10ぐらいのレンジで分布していて、右の軸で表現したい第2変数が0〜100ぐらいのレンジで分布してるとすると、縮尺は1:10になるので、

  1. グラフ領域はまず第1変数にあわせて、y軸が0〜10になるように描く。目盛りは左軸に表示される。
  2. 第1変数をプロットする。
  3. 第2変数を0.1倍して縮め、同じ領域に第2変数のグラフを重ねる。
  4. sec.axisというオプションをつかって、左軸を好きなように変換した軸を右側に設定することができるので、ここで「左軸を10倍する」という変換設定をする。
  5. breaksで目盛りも適当な間隔で設定する。


という手順で、2軸グラフが描かれるわけです。


以下の実行例では、さっき乱数でつくった変数をもっかい使ってるのですが、y1(左軸)y2(右軸)の縮尺を先に計算してscalerという変数に入れてあります。

# データフレームにまとめる
d <- data.frame(Date=x.date, Y1=y1, Y2=y2)

# 各軸の範囲をきめる
# yも最小値〜最大値の形で設定してよいが明示的に与えたい場合が多い気がする
x.lim  <- c(min(d$Date), max(d$Date))
y1.lim <- c(0, 4)
y2.lim <- c(0, 400)


# 左軸と右軸の関係を表すスケーラをつくる
# 各軸の最大最小差の比をとっている
scaler <- (y1.lim[2] - y1.lim[1])/(y2.lim[2] - y2.lim[1])

d %>% 
  ggplot() + 
  geom_bar(aes(x=Date, y=Y1), stat='identity', width=0.7) +
  geom_line(aes(x=Date,y=Y2*scaler, colour = "Y2"), size=1) +
  scale_x_date(date_labels="%m/%d",date_breaks="7 day", 
               limits=x.lim, expand=c(0, 0)) +
  scale_y_continuous(limit=y1.lim, expand = c(0, 0), 
                     sec.axis=sec_axis(trans = ~ ./scaler, 
                        breaks=seq(from=y2.lim[1], to=y2.lim[2], by=50), 
                        name="\nY2\n")) +
  theme(text=element_text(size=12, family="MS Gothic",color="black"),
        axis.text=element_text(size=10, family="MS Gothic",color="black"),
        legend.position="top",
        plot.subtitle=element_text(size=10, color="#666666")) + 
  labs(x="\nDate\n", y="\nY1\n", color = "",
       title='\nggplot2で重ねたグラフ', 
       subtitle='(折れ線を複数追加することもできます)')


f:id:midnightseminar:20200811214800p:plain


sec.axisの中のtransというところには、左軸と右軸の対応関係をformula形式で書くのですが、「.」はデータ全体を表してて、これをスケーラで割るという変換を設定してあります。
x軸が日付なので、scale_x_dateをつかって日付表示の設定をしています。
breaksってところで、右2軸の目盛りを設定しています。水平のグリッドがある場合、左軸の目盛りと右軸の目盛りが噛み合ってたほうが綺麗なので、最初にy1とy2の幅を設定する時に、いい感じの公約数がある値を選ぶのがいいと思います。
themeでフォントを指定してるのは、日本語を文字化けなく表示させるためです。


冒頭に貼った成果物のように、折れ線を2本引きたいときは、geom_line()をもう1個プラスすればいいですね。
グラフのタイトルや軸のタイトルの前後に改行(\n)を入れているのはなんとなく隙間を開けるためです。

*1:だから逆に、重ねる時にスケールが揃ってないことを忘れたりすることがありますねw

Macでスクリーンショットの保存先をショートカットキーで振り分ける

ツイッターで尋ねられて、自分でも気になったので設定してみました。
Macでスクリーンショットを撮るときに、ショートカットキーを工夫して、保存先を振り分ける設定です。
Macのスクリーンショットは、じつはいろいろオプションがあって、ウィンドウ単位で撮影するときに影をつけるかどうかとかも変えられるのですが、私は基本的には、


shift+command+3(全画面を撮影)
shift+command+4(インタラクティブモードで範囲を指定して撮影。space)


の2パターンしか使いません。
で、仮に用途に応じて2つのフォルダに振り分けたいとすれば、計4つのショートカットがあればいいことになります。振り分けるフォルダは、仮に、
"~Desktop/ScreenShot1"
"~Desktop/ScreenShot2"
だとしておきます。


どうやってショートカットをつくるかというと、

  1. 処理自体は、シェルスクリプトで書く
  2. Automatorで、そのシェルスクリプトを呼び出すService(Auromator上ではQuick Actionと呼ばれるのでややこしいが)を作っておく
  3. System Preferences(日本語環境の場合は「環境設定」) > Keybord > Shortcuts > Servicesで、さっきつくったサービスを呼び出すショートカットを設定する*1


という手順になります。
サービスってのは、たとえば↓の画像の場合だとFinderの歯車マークの中に並んでいる、様々な処理のことで、私は正確に理解してないですが「Mac上のいろいろな場面で、ワンクリックで呼び出せる処理」ぐらいに思っています。


f:id:midnightseminar:20200710143305p:plain


さて設定手順ですが、まずAutomatorで、New DocumentをQuick Actionとして立ち上げます。


f:id:midnightseminar:20200710142929p:plain


左のメニューの、Utilitiesから、Run Shell Scriptを、右のワークフローのところにドラッグします。その上にある設定欄は、一行目のinputを設定するところを、no inputにしときます。


f:id:midnightseminar:20200710142952p:plain


次にスクリーンショットを撮って名前を付けて保存する処理をbashで書きますが(上の画像ではもう書いてありますが)、スクリーンショットを撮るコマンドである"screencapture"の細かいオプション設定は下記を参考にしてください。
https://do-zan.com/mac-terminal-screencapture/


ここでは、4パターンのうち、インタラクティブモードで撮影し、ファイルをデスクトップの"ScreenShot1"というフォルダに保存するパターンだけやります。ファイル名は、頭に"sc"をつけて、その後に"20200710140230"(2020年07月10日14時02分30秒)みたいに日付と時刻を並べて、拡張子".png"をつけることにします。
1行目は時刻を取得して格納する処理です。
"-i"ってのが、インタラクティブモードを表すオプション設定。
その後ろは、要するに保存するファイルのパスですね。${current_time}でさっき取得した時刻をくっつけています。

current_time=`date +"%Y%m%d%H%M%S"`
screencapture -i ~/Desktop/ScreenShot1/sc${current_time}.png


これをAutomatorのスクリプトの欄に書いて、たとえば"ScreenShot(Interactive-ScreenShot1)"という名前で保存します。


次に、System Preferences(日本語環境の場合は「システム環境設定」) > Keybord > Shortcuts > Servicesに進むと、一覧の下のほうに、さっきつくった"ScreenShot(Interactive-ScreenShot1)"がありますので、


f:id:midnightseminar:20200710143729p:plain


その右のAdd Shortcutっていうボタンをクリックして、例えば
shift+command+7
と押します。すると、下記のように、ショートカットが登録されます。
(なお、shift+command+6のショートカットは、タッチバーを撮影する機能に割り当てられてて使えませんでした。)


f:id:midnightseminar:20200710143751p:plain


さっきつくったサービス(クイックアクション)を修正したい場合は、
~/Library/Services
という場所にファイルが保存されてるので、それをAutomatorで開いて編集すればいいです。(さっきショートカットを設定した画面から右クリックでもその場所が開けます)


で、さっそく
shift+command+7
を押してみたときに、保存先のフォルダが事前に作成されてないとエラーになります。
あと、場合によっては、スクリーンショットを撮る権限を、Automatorに与えないといけないです。↓のように、システム環境設定の、Security & PrivacyのScreen Recordingの権限を、Automatorに与えます。
これで完了です。


f:id:midnightseminar:20200710143143p:plain


なお新しいサービス(Quick Action)をAutomatorで作成したのに、ショートカット設定画面の一覧に出てこない場合は、システム環境設定のウィンドウをいったん落として開き直すと、出てきたりします。
あと、このやりかただと、まずないとは思いますが1秒以内に2回撮影すると、1回目のは上書きで消えてしまいます。

*1:私はMacの言語設定を英語にしてるのですが、英語の勉強のためでも、カッコつけてるわけでも、ましてや英語のほうが理解しやすいからでもなく、分からないことを調べるときに英語の機能名称とかでググったほうが圧倒的に情報が多いからです。一時期、なぜか英語設定にしてないと無料版がインストールできないソフトとかがあったってのもありますが。

Pythonでよく忘れる、よく間違える書き方

  • Pandasのデータフレームは、=でコピーしようとすると、コピーじゃなくて参照渡しになるので、コピーしたつもりのdfを処理すると元のdfも処理されてしまう。df2 = df1.copy()とするのを忘れないように。
  • Pandasで要素がNaNかどうかを判定させようとする時、要素を指定してその後ろに.isnull()メソッドを使うんじゃなくて、データフレームに.isnull()を実行してそのあとで要素を指定する。
df.isnull().iloc[1,1]
  • リストのappendは、直接編集するメソッドなので、=で代入しなおさなくてよい。(リンク
  • リストをソートする時も、.sort()をつけるだけで、代入はしなおさなくてよい。
  • Pandasから.to_csvするとき、「quoting=csv.QUOTE_NONNUMERIC」という引数をつけると文字列のフィールドは""をつけることができるが、これは事前にimport csvしてないといけない。
  • リストの要素を置換する時、置換対照表を辞書にしておいて内包表記を使うと簡単にかけたりする。
X = ['あ','い','あ','あ','う','え','お','あ']
Y = {'あ':'ア','お':'オ'}
X2 = [Y[x] if x in Y.keys() else x for x in X]
print(X2)
  • listをユニークにするには、NumPyなら.unique()が使えるが、setにしてからlistに戻せばユニークにはなる。順序は崩れるので注意。
unique_list = list(set(original_list))
  • 数値をゼロ埋め(ゼロパディング)した文字列がほしい時は、数字を文字列にしてから、たとえば.zfill(3)とつけるとゼロ埋め3桁の数字(文字列だけど)になる。
  • Rばかりやってると間違えるのだが、Pythonでは代入はオブジェクトのコピーを意味せず、名前が増えるだけとなる。つまり、代入した先の変数名に対して編集すると、元のオブジェクトが編集されてしまう。
x1 = [1,2,3,4,5]
x2 = x1
x2.remove(3)
print(x1)
print(x2)
[1, 2, 4, 5]
[1, 2, 4, 5]
  • andとorを組み合わせると、andが優先(先にくくられる)なので、A and B or C and Dだと(A and B) or (C and D)の意味になる。明示的に優先させたければ()でくくる。
  • nanとnanは==マッチしないので、nanかどうか判定するにはmath.isnan()を使う必要がある。
  • object型の列は、中に複数種類のデータ型が含まれている可能性がある。特に、pandasで文字列型の列を読み込むとobject型になるが、nanを含んでいる場合、そのnanが別の型になってたりとか。


順次追記していきます。

Rでジニ係数を算出

めちゃくちゃ簡単なしょうもない内容ですが、あとで個人的に使うので、ジニ係数を出す簡単な関数をメモしておきます。
ここでは例として、都道府県の人口データのジニ係数を出してみます。
データを小さい順に並べて、都道府県数の累積(これは単順番を表すindexと同じになる)の累積構成比をx軸、人口の累積構成比をy軸にとったのがローレンツ曲線で、そっから三角形とスイカ形の面積比を求めるわけですね。

> # 北海道から沖縄県までの人口(千人単位)が順番に入ったデータ
> print(pop)
 [1]  5286  1263  1241  2316   981  1090  1864  2877  1946
[10]  1952  7330  6255 13822  9177  2246  1050  1143   774
[19]   817  2063  1997  3659  7537  1791  1412  2591  8813
[28]  5484  1339   935   560   680  1898  2817  1370   736
[37]   962  1352   706  5107   819  1341  1757  1144  1081
[46]  1614  1448
> 
> # ジニ係数を取得する関数
> get_gini <- function(v){
+   v.sorted <- v[order(v, decreasing=F)]
+   x <- 1:length(v)
+   y <- cumsum(v.sorted)/sum(v.sorted)
+   AB <- ((length(v))/2)  # 三角形の高さは1なので2乗しない
+   B <- sum(y)
+   gini <- (AB-B)/AB
+   return(gini)
+ }
> 
> # グラフを描く関数
> plot_gini <- function(v, col="black", lty=1, lwd=1){
+   v.sorted <- v[order(v, decreasing=F)]
+   x <- c(0, (1:length(v))/length(v))
+   y <- c(0, cumsum(v.sorted)/sum(v.sorted))
+   plot(x=x, y=y, type="l", col=col, lty=lty, lwd=lwd)
+   par(new=T)
+   plot(x=x, y=x, type="l", col="black", axes=F, ann = F)
+ }
> 
> # ジニ係数を出す
> get_gini(pop)
[1] 0.4415188
> 
> # 図示する
> plot_gini(pop)


f:id:midnightseminar:20200330184908p:plain


斜め45度の線と赤い弧に囲まれた部分の面積が、三角形に占める割合がジニ係数で、1に近づくほど偏りが大きいことになりますね。


ちなみに、1920年、1980年、2018年を比較すると人口ジニ係数は0.248, 0.396, 0.442とどんどん上がっており、図示すると以下のような感じでした。

> # pop1920, pop1980, pop2018はそれぞれの年の都道府県人口
> plot_gini(pop1920, col="#00CC00", lty="dashed")
> par(new=T)
> plot_gini(pop1980, col="#0000CC", lty="dashed")
> par(new=T)
> plot_gini(pop2018, col="#CC0000", lwd=1.5)
> legend(
+   "bottomright",
+   legend=c("1920年", "1980年", "2018年"),
+   lty=c("dashed", "dashed", "solid"),
+   lwd=c(1,1,1.5),
+   col=c("#00CC00", "#0000CC", "#CC0000")
+   )


f:id:midnightseminar:20200330192507p:plain

肺炎死者数の月別データで季節調整を試してみる(RからX-13ARIMA-SEATSを利用)

Rでの季節調整をやってみます。
季節調整は、arima関数のseasonal引数を指定したSARIMAモデルや、stl関数でできてそっちは簡単なのですが(参考1参考2参考3)、一般的にアメリカ商務省センサス局の「X-12-ARIMA」とか「X-13ARIMA-SEATS」とかが有名なので、それを使ってみようと思います。


1965年に開発された「X-11」は移動平均をつかった成分分解法で、これにRegARIMAモデルをつかって事前にデータ補完する機能などを加えたX-12-ARIMAが1996年に開発され、さらに2012年に、X-11とは別系統の技術として利用されてきたARIMAモデルベースの分解方法であるSEATSを選択肢として使えるようにしたX-13ARIMA-SEATSがリリースされた……ということらしい。なのでX-13はX-12が新しくなったというより、X-12とは別の技術も使えるようにしたパッケージという感じのようです。

X-12-ARIMAは、大きく分けて、異常値等を除去するための回帰変数の決定とARIMAモデルによる予測値を推計する「事前調整パート(RegARIMA)」、X-11による季節調整を行う「移動平均パート」、季節調整系列の適切性や安定性に関する診断を行う「事後診断パート」の3つのパートから構成されている。
https://www.stat.go.jp/training/2kenkyu/ihou/71/pdf/2-2-712.pdf

ということなんですが、要するに季節調整のコア部分にはX-11が使われ続けていて、X-12というのはそれに便利機能を付加したものというイメージです。


Rのパッケージだけではできないので、インストール方法からメモしておきます。具体的には、X-13ARIMA-SEATSをインストールして、それをRの{seasonal}パッケージからよびだして使う形になります。{x12}というパッケージからも使えるので追加的に後でやります。
(OSはMacで、X-code、commandline tools、homebrewがインストールされていることを前提とします)


まず、X13-ARIMAのコンパイルにはgfortranというコンパイラが必要らしいですが、これはgccに統合されているとのことなので、gccをインストールすればよい。homebrewで入ります。

$ brew install gcc


次に、センサス局のサイトにいって、X13のソースコードをダウンロードしてきます。
X-13ARIMA-SEATS Download Page (Unix/Linux Version) - X-13ARIMA-SEATS Seasonal Adjustment Program - US Census Bureau
このURLにある、「x13ashtmlsrc_V1.1_B39.tar.gz」という圧縮ファイルをダウンロードして、解凍するとフォルダが出てきますので、それを適当な場所に置きます。私は、
~/Library/R
に置きました。


で、この記事に書いてあるのですが、
Compiling X 13ARIMA SEATS from Source for OS X · christophsax/seasonal Wiki · GitHub
コンパイル前に、makefileの中身を少し書き換えないと、エラーが出ます。さっき解凍したフォルダの中に、
makefile.gf
というファイルがあるので、これをテキストエディタで開いて、293行目の"-static"を消します。


削除前:$(LINKER) -static -o $@ $(OBJS) $(LDMAP) $(LIBS) $(LDFLAGS)
削除後:$(LINKER) -o $@ $(OBJS) $(LDMAP) $(LIBS) $(LDFLAGS)


あとはターミナルからmakeすればいいです。

$ cd ~/Library/R/x13ashtmlsrc_V1.1_B39
$ make -f makefile.gf


コンパイルが終わると、フォルダの中に
x13asHTMLv11b39
という実行ファイルが出てきます。こいつがX13-ARIMAの処理をやってくれるのですが、Rの{seasonal}パッケージは、ディレクトリを指定してやるとその中にある「x13ashtml」というファイルを探すようになっているので、ファイル名を「x13ashtml」に変えておきます。

$ mv x13asHTMLv11b39 x13ashtml


ここまできたら、あとはRのほうから設定を行います。

> library(seasonal)  # パッケージを読み込み
> Sys.setenv(X13_PATH = "~/Library/R/x13ashtmlsrc_V1.1_B39")  # 実行ファイルの置き場を指定
> checkX13()
X-13 installation test:
  - X13_PATH correctly specified
  - binary executable file found
  - command line test run successful
  - command line test produced HTML output
  - seasonal test run successful
Congratulations! 'seasonal' should work fine!


上のように、checkX13()を実行して上記のようなメッセージが出たら、セットが成功していて、X13をRから呼び出せるようになっているということになります。インストールに不具合があった場合、ここでエラーが出ます。


季節調整ができるかどうか、肺炎の月別死者数のデータをつかって試してみます。
データは、以下のURLに載っている「人口動態統計月報(概数)」を取りました。平成17年1月から、死因別の死者数が載っていて、前年分も書いてあるので、平成16年1月から数字が手に入りました。
人口動態調査 結果の概要|厚生労働省


なお、2017年1月から死因の分類変更があって、それまで肺炎に分類されていた人が一部、神経系や認知症系の方に分類されるようになったらしくて、内訳をみて揃えることができなかったので、2017年1月以降は1とするダミー変数も作っておきました。あとで外生変数として使います。


データを読み込んで、とりあえず原系列のグラフ描いてみます。

# CSVの読み込みとts型に変換
# CSVの1列目は年、2列目が肺炎の死者数、4列目に分類変更後ダミーを入力しといた
# (3列目は、ここで使いませんが、自殺者数を入れてある)
d <- read.csv('source/death.csv')
d <- ts(d, start=c(2004,1), end=c(2019,10), frequency = 12)

# 原系列のグラフ
plot(d[,2], xlab="", ylab="肺炎死者数", xaxp=c(2004, 2020, 8))


f:id:midnightseminar:20200318202159p:plain


y軸の目盛りの省略が気持ち悪ければylimを設定すればいいと思います。
これに、季節調整をかけてみます。ひとまずダミー変数を使わずに、死者数の系列をそのまま使います。
X-12も含めて季節調整の原理については『入門 季節調整―基礎知識の理解から「X‐12‐ARIMA」の活用法まで』とかに詳しく書いてあります(ちゃんとよんでないけど)。

まずは推定して結果のサマリをみてみます。

> pneu.x13.1 <- seas(x=d[,2])
> summary(pneu.x13.1)

Call:
seas(x = d[, 2])

Coefficients:
                   Estimate Std. Error z value             Pr(>|z|)    
Constant           0.022586   0.003928   5.751      0.0000000088872 ***
AO2005.Mar         0.217141   0.032563   6.668      0.0000000000259 ***
LS2015.Feb        -0.114776   0.028185  -4.072      0.0000465558118 ***
LS2017.Jan        -0.246416   0.027709  -8.893 < 0.0000000000000002 ***
AR-Nonseasonal-01  0.696167   0.052020  13.383 < 0.0000000000000002 ***
MA-Seasonal-12     0.763174   0.053533  14.256 < 0.0000000000000002 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

SEATS adj.  ARIMA: (1 0 0)(0 1 1)  Obs.: 190  Transform: log
AICc:  2644, BIC:  2665  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 30.71   Shapiro (normality): 0.9937  
Messages generated by X-13:
Warnings:
- At least one visually significant trading day peak has been found in one or more of the estimated
  spectra.


「coefficients」の中にある、「AO2005.Mar」というのは、「Additive Outlier」というタイプの異常値が2005年3月に検出されていて、そのダミー変数の係数が推定されてるわけですね。これはいわゆるスパイク的な変化のことでしょう。「LS2015.Feb」と「LS2017.Jan」は「Level Shift」というタイプの異常値が検出されていて、名前はそのまんまで、水準が階段のように上がったり下がったりする変化です。2017年1月のほうは、分類変更の影響です。なおこれらの他に、増減傾向の変化を表す「TC: Trend Change」というタイプの異常値もありますが、これは今回は検出されませんでした。
「ARIMA: (1 0 0)(0 1 1)」とあるのは、いわゆるARIMA(p,d,q)(P,D,Q)の次数ですね。pはARの次数、dは何階和分過程か(何階階差を取って分析するか)、qはMAの次数を表しています。小文字のほうは時系列方向のARIMAの次数で、大文字のほうは季節方向の次数です。季節方向というのは言い方が難しいですが、年を行、月を列とした表とつくった場合に縦方向に眺めるようなイメージで、たとえば2次のAR(2)だったらその部分はYtがY_t-12、Y_t-24で説明されるということです。
「Transform: log」は、デフォルトで自動変換が選択されてて、ここでは対数変換されたようです。対数変換されているばあい、いわゆる「乗法型季節調整」になり、原系列=トレンド成分*季節成分*不規則成分という分解がされていることになります。


グラフを描いてみます。

# グラフをかく
par(mfrow = c(4,1), oma=c(3,2,2,0), mar=c(0,2,0,2), mgp=c(1,1,0))
plot(pneu.x13.1$series$s10, xlab="", ylab="季節成分", xaxt="n", yaxt="n")  # 季節成分
plot(pneu.x13.1$series$s11, xlab="", ylab="季節調整値", xaxt="n", yaxt="n")  # 季節調整値
plot(pneu.x13.1$series$s12, xlab="", ylab="トレンド成分", xaxt="n", yaxt="n")  # トレンド成分
plot(pneu.x13.1$series$s13, xlab="", ylab="不規則成分", yaxt="n", xaxp=c(2004, 2020, 8))  # 不規則成分
par(default.par)


f:id:midnightseminar:20200318202229p:plain


こんな感じ分解されましたが*1、季節調整値(2段目)やトレンド(3段目)をみると、2017年に急激に下がってますね。さきほどの「分類変更」の影響です。
なお、s10とs12とs13をかけると、原系列に戻せます。また、季節調整値というのは、原系列を季節成分で割ったものになります。(そのへんの話はこの資料をみると分かる。)
ここでは貼り付けませんが、季節成分の抽出結果のテーブル(行が年、列が月)を眺めていると、肺炎死者数はやはり2月や1月が多く温かい時期に減りますが、不思議なことに8月は一貫して7月や9月に比べて増える傾向がある。あと、冬場の肺炎死者数と、4月や9・10月といった春秋の肺炎死者数の格差が、年々小さくなっていました(冬場のウェイトが小さくなって春秋のウェイトが大きくなっている)。


ちなみに、推定結果のオブジェクトをそのままplotに投げると、以下のように、原系列と季節調整値を重ねて、かつ3種類のoutlier(異常値)を図示したものが出てきます。

plot(pneu.x13.1)


f:id:midnightseminar:20200319011509p:plain


予測値と信頼区間もついたグラフを描くならこうします。(描画に凝るならggplot2のほうがいいかも)

### 予測値ありのグラフを描く
# 予測値の取り出し
# 1列目が予測、2列目がlower、3列目がupper
pneu.x13.f1 <- series(pneu.x13.1, c('forecast.forecasts')) 

# 元の時系列のプロット
plot(pneu.x13.1, xlim=c(min(time(d)), max(time(pneu.x13.f1))))

# 信頼区間の塗りつぶし
# 信頼区間下限を左からなぞり、上限を右からなぞる設定でxy座標を作成しpolygonを描画
# time(d)[nrow(d)]とd[nrow(d),2]を入れるのは、実績の最後と
# 予測の最初のあいだに、隙間を開けさせないため
x.ci <- c(time(d)[nrow(d)], time(pneu.x13.f1), rev(time(pneu.x13.f1)), time(d)[nrow(d)])
y.ci <- c(d[nrow(d),2], pneu.x13.f1[,2], rev(pneu.x13.f1[,3]), d[nrow(d),2])
polygon(x.ci, y.ci, col = '#ccccff', border = NA)

# 予測値のプロット
# 隙間を開けさせないために、実績値の最後の値と、予測値の系列をつないで新たなtsオブジェクトにする
pred <- ts(c(d[nrow(d),2], pneu.x13.f1[,1]), start=c(2019,10), end=c(2022,10), frequency=12)
lines(pred, col="#5555ff", lwd=1, lty=1)


f:id:midnightseminar:20200319140341p:plain


次に、2017年1月以降とそれ以前を区別するダミー変数をつかった推定をしてみたいと思ったのですが、なんかエラーが出まくりました。
下記の記事に書いてあるように、外生変数を使う場合、将来予測用に余計なデータが求められて、ややこしいことになるらしい(X13の仕様により)。そこで、x11という引数もあわせて指定すると動くのですが、要するにこれは、成分分解をSEATSではなくX-11で行うということでしょうね。
NOTE not shown · Issue #200 · christophsax/seasonal · GitHub


とりあえずやってみます。
xregという引数に、外生変数を指定します。

> d_sadj.2 <- seas(x=d[,2], xreg=d[,4], 
+                forecast.maxlead = "0", series.span=", 2019.10", x11 = "")
> summary(d_sadj.2)

Call:
seas(x = d[, 2], xreg = d[, 4], forecast.maxlead = "0", series.span = ", 2019.10", 
    x11 = "")

Coefficients:
                   Estimate Std. Error z value             Pr(>|z|)    
xreg              -0.246421   0.027708  -8.894 < 0.0000000000000002 ***
Constant           0.022586   0.003928   5.750      0.0000000089088 ***
AO2005.Mar         0.217145   0.032563   6.668      0.0000000000259 ***
LS2015.Feb        -0.114767   0.028183  -4.072      0.0000465725741 ***
AR-Nonseasonal-01  0.696129   0.052021  13.382 < 0.0000000000000002 ***
MA-Seasonal-12     0.763082   0.053542  14.252 < 0.0000000000000002 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

X11 adj.  ARIMA: (1 0 0)(0 1 1)  Obs.: 190  Transform: log
AICc:  2644, BIC:  2665  QS (no seasonality in final):    0  
Box-Ljung (no autocorr.): 30.71   Shapiro (normality): 0.9937  
Messages generated by X-13:
Notes:
- Unable to test LS2017.Jan due to regression matrix singularity.
- Unable to test LS2017.Jan due to regression matrix singularity.


xregが説明変数として導入したダミー変数で、係数は先ほどのLS2017.Janとほぼ同じ値になってますね。そしてLS2017.Jan は消えています。(一番下のNotesに、重複するので削る旨が書いてあります。)
次にグラフです。

par(mfrow = c(4,1), oma=c(3,2,2,0), mar=c(0,2,0,2), mgp=c(1,1,0))
plot(d_sadj.2$series$d10, xlab="", ylab="季節成分", xaxt="n", yaxt="n")  # 季節成分
plot(d_sadj.2$series$d11, xlab="", ylab="季節調整値", xaxt="n", yaxt="n")  # 季節調整値
plot(d_sadj.2$series$d12, xlab="", ylab="トレンド成分", xaxt="n", yaxt="n")  # トレンド成分
plot(d_sadj.2$series$d13, xlab="", ylab="不規則成分", yaxt="n", xaxp=c(2004, 2020, 8))  # 不規則成分


f:id:midnightseminar:20200318202246p:plain


seriesから取り出す値が「d10」とかになってるんですが、これはX-11で分解(decompositionのdかな?)された結果のデータです。SEATSでの分解結果は、さっきやったように、「s10」とかで取り出される。
トレンド成分をみると、2017年1月の分類変更は取り除かれている感じになっています。「季節調整値」だけだと、単なる季節調整の結果なので、分類変更の影響が出てしまいますが。
なお、上では示してませんが、ダミー変数にかかる係数の推定結果とかも格納されてて、将来予測を出すときはそれも使われることになります。


性能とかは判断できませんが、インストールも面倒だし、X13の仕様を詳しく見ないと使いこなせないので、ほかのもっと基本的なパッケージ・関数を使った方が楽な気はしますねw


ちなみに、↑でインストールしたX13のエンジンは、{x12}というパッケージからも使えます。
外生変数を使わないなら、{x12}のほうがplotの描き方が楽でいいと思いました。

library(x12)
x12path("~/Library/R/x13ashtmlsrc_V1.1_B39/x13ashtml")  # 実行ファイルまで記述

d_x12.1 <- x12(d[,2])

plot(d_x12.1, ylim=c(0, 14200),
     original=T, col_original="black",
     sa=T, col_sa="blue",
     trend=T, col_trend="red", lwd_trend=1.5)


f:id:midnightseminar:20200318224109p:plain


ただ、さっきの分類変更を反映するダミー変数みたいに、回帰の変数を使った推定をしようとするのであれば、設定がいろいろややこしそうです。そもそもX12やX13の機能とシンタックスに関する基礎知識が必要になる模様。
X12やX13の基本が分かっていると、かなりいろいろな要因を季節調整時に考慮できるみたいですが、少し調べたぐらいでは使えそうになかったので、ひとまず放置しました。季節調整を凝らなければならない用事ができたときにまた勉強しようと思います……。


ちなみに自殺者数はこのようになってました。


f:id:midnightseminar:20200318224142p:plain


全体としては減り続けていて、東日本大震災のあった2011年3月に大きく減少し、その後5月に激増しました。

参考

Rの{seasonal}パッケージの仕様
https://cran.r-project.org/web/packages/seasonal/seasonal.pdf

米センサス局のX-13ARIMA-SEATSの仕様
https://www.census.gov/ts/x13as/docX13AS.pdf

X12に関する日本語の説明サイト
http://www.cirje.e.u-tokyo.ac.jp/research/dp/2001/2001cj47.pdf
https://ecitizen.jp/x-12-arima/x-12-arima/

Rの{x12}パッケージの使い方の解説資料(テキストとスライド)
https://www.jstatsoft.org/article/view/v062i02
https://www.r-project.org/conferences/useR-2015/presentations/135.pdf

Rの{x12}パッケージの仕様
https://cran.r-project.org/web/packages/x12/x12.pdf

X-11、X-12-ARIMA、X-13ARIMA-SEATSの関係の解説
http://www.esri.go.jp/jp/archive/snaq/snaq150/snaq150d.pdf

X-11のベースになっている移動平均法による加法型・乗法型の季節調整の原理については、下記資料のp.79あたりからをみると分かりやすい。X-11の問題点*2をX-12がどのように解決するかの説明も書かれてある。
https://www3.boj.or.jp/josa/past_release/chosa199605e.pdf

*1:目盛りの数字が重ならないように設定するのが面倒だったので、y軸の目盛りを消してしまってます。4つのグラフで縦方向のスケールが全然違うのでそこは注意が必要。

*2:前後7年移動平均を取るが末端から7年以内のところではデータがないので後方移動平均に変わっていく点や、妥当性の事後評価ができない点。