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

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

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が別の型になってたりとか。
  • torch.tensorの次元や長さについて。
# 次元数の確認(何重の入れ子になってるか)
x.dim()
x.ndimension()
x.ndim

# 各次元ごとのサイズを確認(外側の次元から順番に要素数が表示される)
x.size()  # numpy.arrayだとこれでバラしたときの要素数が返る?ので注意
x.shape

# 要素数(次元をバラしたときの要素数)
x.numel()
x.nelement()
  • 文字列の出力を、下に並べていくんじゃなくて上書きさせたいときは、キャリッジ・リターンで先頭に戻す動きをすればよい。(前の出力のほうが長くてケツのほうが残ったりするのを避けるには、スペースで埋める。)
for i in range(100):
    print(f'{i}週目です', end='\r')
    time.sleep(5)
  • データフレームを複数の条件でソートするとき、Rと違って、高階層側(外側)の条件を先に書く。以下のようにすると、全体としてはmatch_id順に並んでいて、同じmatch_idの中でplayer名順に並ぶようになる。Rとかエクセルだと内側から先にソートするので注意。
d2 = d2.sort_values(by=['match_id','player'], ascending=True)
  • pandasのデータフレームをソートしたあと、reset_index(drop=True)しておかないと、順番は代わってるのにインデックスのラベルはもとのままになっているので、あとでおかしなことになる。
d2 = d2.sort_values(by=['match_id','player'], ascending=True)
d2 = d2.reset_index(drop=True)

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年以内のところではデータがないので後方移動平均に変わっていく点や、妥当性の事後評価ができない点。

四半期GDP成長率の「年率換算」はどれぐらいブレるのか?

GDPは「四半期速報値」というのが作成されており、前の四半期との比(たとえば2019年7-9期と10-12月期の比)を4乗したものが、「年率換算」の成長率として報道などで使われている。最新の速報では、2019年10-12月期の実質GDPが、年率換算でマイナス7.1%と大幅な減少を記録したと報道されていた。


この、四半期成長率を4乗した値というのは、トレンドが安定していればもちろん年率の成長に一致するが、短期的な変動が4乗されるわけなので、けっこう大きく上下する。

以下は、ちょっとわかりにくいが、実質GDPの
-年間(暦年)の成長率[黒い線]
-四半期成長率を「年率換算」したもの(前期比の4乗)[赤い線
-四半期GDPを4倍したものと、前年の年間GDPの比[青い線
をグラフにしたものだ。四半期はいずれも季節調整済み。


f:id:midnightseminar:20200309150830p:plain


データは、暦年GDPについては平成30年度の国民経済計算年次推計、四半期速報値は2019年10-12月期の一次速報値(実質・季節調整済み)(年率換算マイナス6.3%と報道されたやつ)を用いた。


年成長率(黒線)と、四半期成長率の年率換算値(赤線)は、けっこう乖離があるということがわかる。ちなみに、「年成長率」と「四半期成長率の年率換算」の相関係数をとると0.33ぐらいになる。

最新の速報では2019年10-12月期の年率換算成長率がマイナス7.1%になったということだが*1リーマンショックの影響が出た2009年1-3月期は年率換算で-17.7%、今回と同じく消費増税があった2014年4-6月期は年率換算-7.4%だったので、今回はまあ、それらよりはマシということになる(それでも大きな落ち込みだが)。次の2020年1-3月はコロナの影響があるので、リーマン並かそれ以上いくのかも?


一方、「四半期GDPを4倍したもの」と「前の年のGDP」の比をとったのが上のグラフでいうと青い線なのだが(「前年同期比」ではない点に注意)、これは年成長率とあまり大きくは乖離していない。年成長率との相関係数も、0.93ぐらいになる。


もちろん、短期的の成長を瞬間風速的にみたい場合はもちろんあるだろうから、「四半期成長率の年率換算」に意味がないわけではない。というか、「年成長率」ですら、人によっては「短期的すぎる」指標かもしれないので、まあ要するにみる人の関心次第としか言いようがない。また、上のグラフでいう青い線のほうが優れた性質のデータだ、というわけでもない。赤と青は意味合いが異なるわけで、単に「年成長率との近さ」を知りたいなら、青い線のほうをみたほうが近くなるという話である。


重要なのは、四半期成長率の年率換算というのは、あくまで「その3ヶ月と前の3ヶ月のあいだの短期的な変動」を誇張したものなので、「年率換算」したからといって、「年スケールでのトレンド」を見通そうとしているわけでは全くないということ。まぁ当たり前だが。

上のグラフで確認したかったのは、「四半期スパンの変動と年スパンの変動はどれぐらい乖離しているのか」であって、四半期GDPをみることが役に立つか立たないかということではない。

ただ、年成長率とはあまり関係ないものを見てるのであり、4乗したところで年成長率に近くなってるわけでもないのだから、換算しなくていい気はするが。四半期の速報値から年スパンでの変化を推し量りたいのであれば、上のように前の年の年間GDPと比較するか、過去4四半期の平均をその前の4四半期の平均と比較するかすればいいんじゃないだろうか。

 

会社員時代、営業データを週次・月次・年次で報告する仕事をしてたのだが、その経験からすると、上のグラフの赤い線(四半期GDPの年率換算)のように大きくジグザグする指標は、直接モニタリングするのにはあまり向いていない可能性もある。時系列データをどういう単位でみるべきかというのは、考え方はいろいろあり、もちろん短期的なショックを検出することも大事ではある。しかし、大きく下がったあとに反動で大きく上がるような数字は、天候とか平日日数*2みたいな特殊要因に大きく左右されてるわけで、解釈がめんどくさい。というか「短期的な解釈」だけをして「トレンドとして解釈」してしまわないような注意がいる。トレンドはトレンドで別途確認するという心がけが必要だろう。

*1:上のグラフは、最新の速報の前の、「マイナス6.3%」になってる1次速報データで作っている

*2:これは四半期GDP統計では考慮済みらしいが

Macのパーティションがおかしくなった場合の対処

Mac OS(最近はmacOSというらしいが)のディスクユーティリティでパーティションを変更しようとするとき、よくわからない現象に直面することがたまによくあります。
下記の記事のように、一部のボリュームが消してくても消せない(マイナスボタンがない)状態になったりします。
mac GUIで結合できないパーティションを結合する – MY ROBOTICS


今回、外付けHDDのパーティションをいじろうとしたらそれに陥ったので、ターミナルからdiskutilコマンドで操作してみたのですが、それでも変なことがおきました。
具体的には、

$ diskutil list

/dev/disk2 (external, physical):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:      GUID_partition_scheme                        *4.0 TB     disk2
   1:                        EFI EFI                     209.7 MB   disk2s1
   2:                  Apple_HFS TM_MBP2016_Portable     1.0 TB     disk2s2
   3:                  Apple_HFS Untitled                699.9 GB   disk2s4
   4:                  Apple_HFS Strage_main             2.3 TB     disk2s5


というふうに、確認すると(diskutil listの出力結果はexternalのところだけ貼ってます)、disk2という外付けHDDの中にEFI(システム側で必要なやつ?)を入れて4つのパートがあって、足したら4TBになります。これのうち、disk2s5だけ残してdisk2s2とdisk2s4を削除し、結合して(EFIを含めて)4TB使えるようにしたかったわけです。要するにパーティションをやめて全部結合するということ。


そこで、

$ diskutil eraseVolume "Free Space" not disk2s2
$ diskutil eraseVolume "Free Space" not disk2s4


というふうにボリュームの消去をしてみたのですが、その結果どうなったかというと、

/dev/disk2 (external, physical):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:      GUID_partition_scheme                        *4.0 TB     disk2
   1:                        EFI EFI                     209.7 MB   disk2s1
   2:                  Apple_HFS Strage_main             2.3 TB     disk2s2


こういうふうに、不要なボリュームは消えてくれました。ディスク全体としては4.0TBあって、残存するボリュームの合計が2.3TB+209.7MBであり、残りがFree Spaceになったのだと思います。
で、現在2.3TBになっているdisk2s2をリサイズしてマックスまで使えるようにしたかったのですが、ディスクユーティリティで見てみると、


f:id:midnightseminar:20200229171848p:plain


こんなふうに、Free Spaceは表示されなくて、全体で2.3TBしかないよということになってしまいました。
そこで、ネットでみた説明を参考にdiskutil resizeVolumeというコマンドをつかって、リサイズ後のサイズを「100%」にしてやってみたのですが、

$ diskutil resizeVolume disk2s2 100%
Resizing to full size (fit to fill)
Started partitioning on disk2s2 Strage_main
Error: -69743: The new size must be different than the existing size


こういう結果になりました。100%だともとのサイズと同一という意味であって、ディスクをマックスまで使うという意味じゃなかったんですねw
なので、具体的に数値で指定してやってみたのですが、たとえば3.9TBに設定したとしても、

$ diskutil resizeVolume /dev/disk2s2 3.9TB
Resizing to 3900000000000 bytes
Started partitioning on disk2s2 Strage_main
Error: -69519: The target disk is too small for this operation, or a gap is required in your partition map which is missing or too small, which is often caused by an attempt to grow a partition beyond the beginning of another partition or beyond the end of partition map usable space


となり、容量が足りませんと出ます。3.9TBにEFIの分を足しても超過するわけではないはずなのですが・・・。
結局、もう調べるのがめんどくさくなったので、
Macのディスクユーティリティーコマンドで外部ストレージをフォーマット - ytooyamaのブログ
を参考に、論理ボリュームではなく物理ディスク全体を消去することにしました。データが全部きえますが、別のディスクにコピーしてあるので後で戻します。

$ diskutil eraseDisk JHFS+ Untitled /dev/disk2


そしたらまぁ、


f:id:midnightseminar:20200229172118p:plain


こんなふうに、4TB使える状態にはなりました(画像はボリューム名をUntitledからStrage_mainに戻したあとのものです)。
Free Spaceと結合して論理ボリュームのサイズをでかくする方法は、後でわかったら追記します。