二次元の散布図で一週間ぐらい盛り上がっているが、玉石混合の意見が流れていって教育上（誰の？）よろしくない。
相関係数が低い非線形の関係の2変数の散布図に、直線をフィットするなという非難が理工系の人々から出ているのだが、統計学的には問題はない。非線形モデルが単調関数であれば一次近似として線形回帰は機能する*1し、直線の傾きの係数の標準誤差が計算されて表示されているからだ。
断り無く除外された観測値がありチェリーピッキングだいう非難もあったが、異常値の除去はよく見られる方法である*2し、同様のデータセットをつくってロバスト回帰で分析しなおしたところ概ね話は再現できた*3ので、研究不正とは言い難い。
相関なのに因果であるかのように説明した、雑な散布図をもとに政策提言したと言う批判は、どうも財務省の紹介文が悪く、それで誤認した人の話が広がったのが強い理由のようだ*4。微妙なキャプションもあるのだが、そちらは微妙。
非難の多くは誤解もしくは「お気持ち」だ。もっとも話題の散布図が無問題（であった）かというと、そうではない。訂正はされたが、データセットの整備と分析対象の抽出方針が雑である蓋然性が高い。あまり言葉にされない（データクリーニングガガガガ…という悲鳴は聞こえる）作業なのだが、かなりのやらかしはここで入る。
1. データセットの整備
自分で個票を集めるのは大仕事で、学術論文にしないような話では、政府や中央銀行、国際機関のデータをダウンロードしてきて編集する事が多いのだが、大概、データ分析に使いやすい並びにはなっていない。話題の散布図の元データは、
と言うようなExcelファイルのシートで、年ごとに横方向に伸びていくものだ。ヒトが目視することを考えている。
しかし、Rで統計解析をする場合は縦横逆の方が扱いやすいし、複数のファイルの複数のシートをつなぎあわせて、以下のようにする必要がある。
| CODE | Country | YEAR.PSF | PSF | YEAR.TFR | TFR | CLASS | YoA |
|---|---|---|---|---|---|---|---|
| AUS | Australia | 2015 | 2.69 | 2017 | 1.74 | H | 1971 |
| AUT | Austria | 2015 | 3.18 | 2017 | 1.52 | H | 1961 |
| BEL | Belgium | 2015 | 3.2 | 2017 | 1.65 | H | 1961 |
| CAN | Canada | 2015 | 1.77 | 2017 | 1.55 | H | 1961 |
| CHE | Switzerland | 2015 | 2.11 | 2017 | 1.52 | H | 1961 |
| CHL | Chile | 2015 | 1.72 | 2017 | 1.56 | H | 2010 |
| COL | Colombia | 2015 | 2.03 | 2017 | 1.72 | UM | 2020 |
| CRI | Costa Rica | 2015 | 1.04 | 2017 | 1.71 | UM | 2021 |
| CZE | Czechia | 2015 | 2 | 2017 | 1.69 | H | 1995 |
| DEU | Germany | 2015 | 3.03 | 2017 | 1.57 | H | 1961 |
| DNK | Denmark | 2015 | 3.53 | 2017 | 1.75 | H | 1961 |
| ESP | Spain | 2015 | 1.44 | 2017 | 1.31 | H | 1961 |
| EST | Estonia | 2015 | 2.94 | 2017 | 1.59 | H | 2010 |
| FIN | Finland | 2015 | 3.1 | 2017 | 1.49 | H | 1969 |
| FRA | France | 2015 | 3.7 | 2017 | 1.89 | H | 1961 |
| GBR | United Kingdom | 2015 | 3.58 | 2017 | 1.74 | H | 1961 |
| GRC | Greece | 2015 | 1.29 | 2017 | 1.35 | H | 1961 |
| HUN | Hungary | 2015 | 3.38 | 2017 | 1.51 | H | 1996 |
| IRL | Ireland | 2015 | 1.98 | 2017 | 1.78 | H | 1961 |
| ISL | Iceland | 2015 | 3.28 | 2017 | 1.76 | H | 1961 |
| ISR | Israel | 2015 | 2.11 | 2017 | 3.11 | H | 2010 |
| ITA | Italy | 2015 | 1.92 | 2017 | 1.34 | H | 1962 |
| JPN | Japan | 2015 | 1.61 | 2017 | 1.43 | H | 1964 |
| KOR | Korea, Rep. | 2015 | 1.28 | 2017 | 1.05 | H | 1996 |
| LTU | Lithuania | 2015 | NA | 2017 | 1.57 | H | 2018 |
| LUX | Luxembourg | 2015 | 3.17 | 2017 | 1.39 | H | 1961 |
| LVA | Latvia | 2015 | 2.27 | 2017 | 1.69 | H | 2016 |
| MEX | Mexico | 2015 | 0.99 | 2017 | 2.04 | UM | 1994 |
| NLD | Netherlands | 2015 | 1.75 | 2017 | 1.62 | H | 1961 |
| NOR | Norway | 2015 | 3.37 | 2017 | 1.62 | H | 1961 |
| NZL | New Zealand | 2015 | 2.63 | 2017 | 1.81 | H | 1973 |
| POL | Poland | 2015 | 1.92 | 2017 | 1.45 | H | 1996 |
| PRT | Portugal | 2015 | 1.51 | 2017 | 1.38 | H | 1961 |
| SVK | Slovak Republic | 2015 | 1.94 | 2017 | 1.52 | H | 2000 |
| SVN | Slovenia | 2015 | 1.8 | 2017 | 1.62 | H | 2010 |
| SWE | Sweden | 2015 | 3.48 | 2017 | 1.78 | H | 1961 |
| TUR | Turkiye | 2015 | 0.38 | 2017 | 2.08 | UM | 1961 |
| USA | United States | 2015 | 1.11 | 2017 | 1.77 | H | 1961 |
4ファイルの内容をつなぎあわせている。なお、PSFが家族関係公的支出（Total）、TFRは合計特殊出生率だ。余談だが、Rで上のテーブルをデータフレーム df01 に入れて、summary(lm(TFR ~ PSF + as.factor(CLASS), df01, subset = CCODE != "ISR")) としたら、中位国を入れても話が変わらない事が確認できる。皆さんが目視で推測している相関係数も出てくる。
Rのコードで機械的に上のように並び替えることはでき、実際にほとんどそのように行った。話題の散布図は2015年と2017年のデータを使っているが、年を変えて試すこともしたいし、差分モデルにするには4項目追加しなければならない。2011年から2020年まで連続してOECDに加盟していてClassがHというような条件もつける。
生成AIの発達で状況はだいぶ変わっているであろうが、コーディングに慣れていない社会科学分野の研究者は多い。コピペ編集の繰り返しは思わぬミスを生むことがある。
もっともプログラムで処理しやすいモノではない。同じ調査のはずなのに項目が変更になることも多い。大学などの進学率の確認で、文科省を呪った人は多いであろう。こっそり項目の定義が変わることすらある。データの跳ねに頭を悩まし、官公庁に理由を訪ねにいったところ…と言う人もいるはずだ。
泥縄になりがちなので、作業記録はつけておく方がよい。ダウンロード元が消えたり、日銀がデータ項目の管理番号を変えたりしてくるが、記録がないと具体的にどこから引っ張ってきたデータか分からなくなる。
2. データセットの正確性チェック
手間暇だけで済むかと言うと、しっかり罠が仕込まれていることが多い。カントリーコードがなく表記にブレがある（e.g. Korea, Rep.とKorea）、欠損値（NA）の表記にブレがあると言うような事が頻繁に起きる。機械的に実行すると、この辺りのブレで問題が入り込む。
コードが生成するテーブルの行数をチェックするのは必須で、幾つかランダムに個体を選んで元データと数字があうか調べることも必要だ。欠損値（N/A）はコーディングミスが由来なことも多いが、原データからしてと言うこともある。なぜか数字の列が文字列扱いになり、数値に強制変換をかけても2011年のデータにNAが出るので原データを見たら、*3.4と言うようにデータに脚注を入れてあったことなどがある。毎月勤労統計のことだが。
3. データの抽出方針に理由をつける
観測値は概ね個体（e.g. 国, 地方自治体, 個人）と時点（e.g. 年, 年度, 月）で決定されるが、それぞれどうするか決める必要があり、文に含めるかは別として、明確な説明をする準備が要る。
個体
- 分析対象にする基準を正しく言葉にする
- 国際比較をすると、国ごとのデータの精度がマチマチで、やはり開発途上国の数字は信用ならない面がある。そもそも数字の公表がないことも。先進国、例えばOECD加盟国に限ることに合理性はある。しかし、「OECD加盟国を分析対照にする」という方針だけでは、何年の時点の加盟国なのか曖昧だ。実際、話題の散布図では分析年の2015年にはまだ加盟していないラトビアとリトアニアが入っている問題が生じている*5。
- 異常値として除去する基準を正しく言葉にする
- 話題の散布図で問題になったイスラエルのような異常値の取り扱いは難しい。宗教的・文化的に…と言うのは容易いが、具体的にそれが何なのか激詰めされる恐れがある。イスラエルの宗教的・文化的特徴を述べた上で、それが変数の関係にどう影響するのか説明する準備が必要だ。分析手法が複雑になっても、正体不明の要因をコントロールする方が楽かも知れない。
- 逆に中所得国のようなカテゴリー全体を排除するのを正当化するのは容易い。多変量解析では、大概、中所得国ダミーを入れるだけで、高所得国だけのサブサンプルと同様の結果を得ることができる。所得による影響があるが、それを入れると三次元になってしまい図示できないので弾いたが、多変量解析では含めることができると説明できる。
時点
- 分析時点を選んだ理由を正しく言葉にする
- 社会データは様々な要因に影響される。長期平均をとるような操作がない場合は、戦争や疫病などの他の強力な要因が無い時点を選ぶ必要がある。
- 一方で、政策論では新しいデータの方が好まれる。時代の変化もあるから、昔の法則が現在も当てはまっているのか分からないからだ。もっとも西南戦争 (1877年)後のインフレーションの話を古すぎると言って拒絶する人は見たことはない。
- 話題の散布図では、それが収録された著作の発売年2019年で入手できる中で最新のデータを用いた蓋然性が高い。国際機関がまとめるデータは、何年も遅れて出てくるからだ。2026年の現在では、もっと新しい数字を見たくなるが、2019年末にはじまったコロナ禍の影響が高い一方、（この記事を書いている時点では）家族関係公的支出は2021年の数字が最新だ。
4. 横軸と縦軸に何をとるのか考える
横軸に説明変数、縦軸に従属変数の値をとるのは自明だが、変数に選択肢がある場合、多変量の関係である場合は考えることが増える。また、どの変数をどう使ったかは明記しなければならない。
話題の散布図では横軸の従属変数となった家族関係社会支出は、Total、Cash、Service、Tax-breaks for familiesの4項目がある。大項目のTotalを使っている蓋然性が高いが、どうも説明されていない。現金給付よりも現物支給というような議論がされていたわけで、細目（もしくは複数の細目の合計）を使うという選択肢もあった。
予測値は一般線形回帰（OLS）のものが通例だが、小技として異常値が多い場合はランク基準回帰などのロバスト推定もありえる。
y = ꞵ₀ + ꞵ₁x + ꞵ₂z + ϵのような多変量回帰の結果をプロットしたい場合は、縦軸にy - ꞵ₂zをとり、横軸にꞵ₀ + ꞵ₁x + ϵをとるような事もアリだが、調整されたyなる概念を導入しないといけないため、場所によっては聴衆がついてこない。
3次元
バブルプロットなどを使うことで、縦軸横軸にそれぞれ異なる説明変数を割り当て、三次元化する選択肢もあるが、これも聴衆受けは同様。丸が大きくなると、他の点と被って見づらくなる。
従属変数が質的データのとき、予測値に応じて（等高線プロットを使って）境界描き、観測値に応じてマーカーの形状や色を変えて点を描き込むのは、ぼちぼち理解してもらえるようだ*6。
5. 図に添える文字情報を考える
散布図に限らず、ソースや編集方法、加工方法などの注記を入れないといけない。役所の図（内閣府の生産関数を推定しないと出てこないGDPギャップ）で、プロットよりも文字情報の方が多いぐらいのものを見たことがあるのだが、文句がつけられたときの抗弁の準備は念入りに行った方がよい。一年ぐらいしたら何をやっていたのか自分でも忘れる。
タイトルは難しい。話題の散布図はあちこちで使い回されていたのだが、使っている場所ごとにタイトルが異なっていた。相手にあわせたか何かだと思うが、誤解を招くようなタイトルに変化してしまっていた。2019年の著作では「家族関係社会支出が高い国ほど出生率も高い」というキャプションだったが、2020年のプレゼン資料では「家族関係社会支出と出生率」に転じた後、2023年のプレゼン資料では「少子化対策が充実した国ほど高出生率」に（Twitter）。家族関係社会支出の中に少子化対策も含まれるが、高齢者支援なども入っているし、コロナ禍支援の現金給付も含まれるので、2021年のデータがアレだと言うのに…(´･ω･`)ﾁﾘﾁﾘﾁﾘ*7
6. 欠損データ処理
欠損データがある行は除去しても良いのだが、除去される行が半分近くになってくると、よくない推定量になる可能性が増す。欠損が生じる確率と、変数の値が相関しているかも知れない。線形回帰であったら今やそう手間暇がかからないので、欠損値補定*8を行う方が望ましい。2変数の散布図のときは、プロットしない変数を、欠損値補定のための補助変数として使いたいところではあるが。
まとめ
理工系ラボの統計解析では、あまり気にしなくて良いことが多数あると思う。すべて大した話ではないのだが、逆にサボりやすい面もあるはずだ。ラボでの実験時には安全面もあって防塵メガネとマスクとヘアキャップと白衣と手袋をしろと指導してコンタミネーションを防ぐわけだが*9、データセット作成の事務作業で死ぬことはない。
なお、二次元プロットと言えども線を引くということは計量分析になっており、小難しい問題に直面することもある。時系列データのプロットで、各点に信頼区間をつけて問題を指摘されている人がいた。状態空間モデルを使うべきだと言うような話であったと思うが、伝わってはいなかった。
*1自然冷却の実験データは非線形モデルに従うことが知られており、ニュートンの冷却法則は線形回帰よりも良い説明力を持つ（ある学部生の実験データで0.9693）が、線形回帰をしても増減の方向を誤るようなことはなく、そこそこの予測精度（相関係数0.7693）を持つ。非線形モデルは推定が頑強ではないことが多く、現象をよく説明するモデルが単調関数であることぐらいしか仮定できない場合は、線形回帰モデルを用いることに合理性がある。
*2プロットすればイスラエルが他と大きく異なるのは一目瞭然だ。異常値も含めて全体の傾向を言うべきではないかと言う批判はありえるが、よくある操作。不正とは言い難い。もっとも異常値の発生を含めて真のモデルという捉え方もできるので、異常値があるときはロバスト推定を使う方をオススメする。
*3関連記事：あるにはあるよ、家族関係社会支出と出生率の相関
*4関連記事：皆さんがあの散布図で東大教授が因果を主張したと受け取ってしまったのは、財務省が悪い
*5この2カ国の有無は推定結果に大きな影響はないので、意図的なものではない。
*6やりたい場合は「Rの等高線プロットで花びらを描く」などを参照されたし…と生成AIの時代であった。
*7チリはコロナ禍で大量の給付金を出し、家族関係社会支出がGDP比で0.72%（2019年）から、7.34%（2021年）に急増した。当時は暴動があったり深刻な状況であった（AFPBB News）。他の国も給付金を出しており、2020年頃に家族関係社会支出と出生率の関係は大きく変化した。
*9安全ピペッターなしのホールピペットで毒物を扱っていたという武勇伝は語らなくてよい。
0 コメント:
コメントを投稿