データ分析メモと北欧生活

旧Untitled Note. データ分析、計量経済・統計とR、水産管理、英語勉強、海外生活などについて備忘録や自分の勉強のOutputの場所として

MENU

【ノルウェー】SbankenでBankIDが発行されない問題(銀行口座開設)

昔書いたこの記事の続きを書き忘れていたことに気がついた。

keita43a.hatenablog.com

問題：本人認証できないので口座開けない。

Sbankenで銀行口座を開こうとすると、BankIDが取得できない。（その後、口座すら開けなくなった）

ノルウェーでの公的なサービスを利用するためにはセキュリティのためにいくつかのIDが存在する。これは個人識別ID(Fødselsnummer)と紐付いているが、別物だ。

自分の税金の状態にアクセスする、ノルウェーの病院に予約を取ったりメッセージを送るHelsenorgeのアカウントにアクセスする、などのためにはセキュリティレベルの高いIDが必要で、その一つがBankIDである。BankIDは銀行口座を開けると取得でき、BankIDは一人に一つなので、いくつも口座を開いても複数のBankIDを取得できるわけではない。逆に、BankIDを持っていると別の銀行で口座を開くときに手続きが早い。

その他にもいろんなシーンで使うことがあるので、ノルウェーに住むなら必須のIDである。IDと言って物理的にカードなどが発行されるわけではなく、個人識別IDに紐付いたパスワードと銀行のワンタイムパスワード発行機で発行するパスワードを使うが、BankID取得後はスマホに登録して使えるBankID på Mobileがあるので、ワンタイムパスワード機を携帯する必要もなくなる。

店舗のないSbankenの問題。

流行りのネットバンクで手続きが早いと聞いていたのでSbankenで開いたのだが、問題は店舗のないことだった。実店舗がないので、本人の認証を郵便局(Posten)で行うというシステムになっている。アカウントを開くときに指示されたとおり、最寄りのPostenでIDをスキャンするのだが、当然IDとしては日本のパスポートしか無い。しかし、郵便局員がパスポートを何度やっても認識されないと言われる。

日本のパスポートはRFID入っているのに、と思ったけど、どうやらスキャン機自体が、ノルウェーで2010年4月以降に発行されたパスポートでないと認識されないらしい。

私が開いたとき(2018年秋）には、スキャン機が入れ替わった直後だったらしくて、郵便局員いわく”古い方法”で、物理的にパスポートのコピーを取って送付することで銀行口座自体は開けたが、それでもBankIDを取得するためにはスキャンが必要らしい。
また私の妻が後に開こうとした2020年の4月ごろは、その方法でも口座自体開くことができなくなっていた。

解決法：結局、実店舗のある銀行で開いた

ではどうするか？と、聞いたところ、公証人役場（ノルウェーでは裁判所？）に行き、パスポートが本人IDをであることを証明してもらう書類を作ってもらい、それを送れ、と言われた。

以下の記事でも、同様に苦労して公的書類を作った学生のエピソードが書かれている。
Exchange student struggles to open bank account - The Oslo Desk

すごくめんどくさそうだったので、結局DNBで作ってしまうことにした。
DNBならオンラインで申請したあと、実店舗でパスポートで認証できる。受け付けてからの時間はかかるが面倒が少ない。(BankIDを持っていれば手続きが早いらしい。BankIDを取得するための手続きなので仕方がない。。。)
留学生などはSpare Bank 1で作ってる人も多い印象。

他にどの銀行がいいのか？というのはこちらの記事も参考にご覧ください。

keita43a.hatenablog.com

【mac】Homebrewをインストールしたのに「command not found」が出る。

Homebrewをインストールしたのに、brewと打ってもコマンドが見つかりませんになったときの対処のメモ。

検索すると、「Homebrewをインストールしてないから」という回答が日本語でも英語でも出てくるのだが、そうじゃねえんだよと。

lanchesters.site


問題は、Homebrewをコマンドラインでインストールできたはずなのに（brew helpと打ってみよう！と言っているのに）、いざ打ってみるとcommand not foundと出てくるということである。（スクリーンショット）

f:id:keita43a:20200617211602p:plain


さて、注目すべきはWarning: /usr/local/bin is not in your PATH.という部分であった。
つまりパスが通ってなかったので、それを通してやればよかった。

keita$ export PATH=/usr/local/bin:$PATH

これでbrewが通るようになった。

参考：
github.com

【R】geom_sfで描画エラーが出たけどsfアプデしてないだけだった

R R 備忘録 データ分析

簡単な備忘録。tidyverseやR本体などアップデートしたら、付随する他のパッケージもアップデートしましょう。

以前は普通に走ってたポリゴンファイルを読みこんでggplotで描画するコードにエラーが出たので調べたら、
tibble (というかtidyverse）はアップデートしたのにsfはアップデートされてなかったのでしましょう、というオチだった。


github.com

エラー

Error: All columns in a tibble must be vectors.
x Column `geometry` is a `sfc_GEOMETRY/sfc` object.
Run `rlang::last_error()` to see where the error occurred.

コード

jp_sh = readRDS("gadm36_JPN_1_sf.rds")

jp_sh2 = st_simplify(jp_sh, preserveTopology = TRUE, dTolerance = 0.01) 

ggplot(jp_sh2) + 
  geom_sf()

sfパッケージをアプデしたら普通にできました。

f:id:keita43a:20200607185220p:plain

他の可能性

他には、開発バージョンのtibble使ってたら同じエラー出たケースとかあるので、tibbleとsfの整合性が問題になる部分があるみたいです。

github.com

【R】新しい回帰分析表のパッケージ {modelsummary}

R データ分析 計量経済学

f:id:keita43a:20201209204033j:plain

Rで計量経済や統計分析やる時に、結果をきちんと理解しながらモデル作りたいし、できたモデルの結果を書き出すのも間違いなくやりたいですよね。

新しい回帰分析表のパッケージを発見したので、ざっと試してみました。

Rの回帰分析のパッケージ

Rでの回帰分析の表を書き出すのには {stargazer}というパッケージがあって、かなり柔軟にかつ自動的にhtmlやLatexに書き出しができるのですが、作者の都合か、更新が止まってしまってます。更新が止まってしまうと、新しい分析パッケージの結果とかを受け付けられなくなるんですよね。特に計量経済における回帰分析で役立つ{estimatr}の結果を受け付けないのがつらいです。

keita43a.hatenablog.com


他のパッケージで書き出すことも可能です。例えば{huxtable}パッケージのhuxreg関数は、かんたんに書き出せるし、Rのコンソールに見やすい結果を出力してくれるので良いです。しかし、まだ発展途上ということで、{stargazer}ほどの柔軟性はないのと、Latexで出力する際に表の位置が真ん中にいかないという問題がありまして、まだ実用的ではないかなという個人的な印象です。

以前なんとかestimatrの出力結果を固定効果付きでhuxregで書き出すという記事も書きました。

keita43a.hatenablog.com

今のところは作業中はhuxregで結果見ながら作業して、最終的に出力するときはいろんなワークアラウンドを駆使してstargazerできれいな表を書き出すという形をとっています。

{modelsummary}パッケージ

そんな時に、以下のTweetを見つけました。

新しい統計モデルのサマリー表ですって？？
しかも、スレッドの次のTweetを見てみると、なかなかきれいな表が出力されています。


パッケージの説明はこちら。

vincentarelbundock.github.io


さっそくインストール

{modelsummary}がCRANに上がっているというツイートだったので、CRANからダウンロードしましょう。
dplyrの最新版に依存してるみたいなので、tidyverseも一緒に更新しておきましょう。

install.package("modelsummary")
install.package("tidyverse")

シンプルに回帰分析の出力

とりあえずmtcarsデータ使って、回帰分析を行って、結果を出力してみます。
結果はmsummary()関数を使います。

reg <- lm(mpg ~ cyl+disp, mtcars)

msummary(reg)

この関数を使うと、コンソールではなくViewerにきれいなテーブルが出力されます。

f:id:keita43a:20200529203525p:plain


２つ以上のモデルを同時に表示させる場合には、summary()のようにオブジェクトをそのまま並べるのではなく、リストにします。
リストの要素指定時にモデル名をつけることができて、これが列名になるようです。

regs <- list()
regs[['Model A']] <- lm(mpg ~ cyl+disp, mtcars)
regs[['Model B']] <- lm(mpg ~ cyl+disp+hp, mtcars)

msummary(regs)

f:id:keita43a:20200529203629p:plain

出力フォーマット

コンソールに表示

オプションでmarkdownと指定すれば、huxregのようにコンソールに表示することも可能です。

msummary(regs, 'markdown')

|            |Model A |Model B |
|:-----------|:-------|:-------|
|(Intercept) |34.661  |34.185  |
|            |(2.547) |(2.591) |
|cyl         |-1.587  |-1.227  |
|            |(0.712) |(0.797) |
|disp        |-0.021  |-0.019  |
|            |(0.010) |(0.010) |
|hp          |        |-0.015  |
|            |        |(0.015) |
|Num.Obs.    |32      |32      |
|R2          |0.760   |0.768   |
|Adj.R2      |0.743   |0.743   |
|AIC         |167.1   |168.0   |
|BIC         |173.0   |175.3   |
|Log.Lik.    |-79.573 |-79.009 |

出力フォーマットは以下のものに対応しているようです。

フォーマット オプション
HTML 'html'
LaTeX 'latex'
Text / Markdown / ASCII 'markdown'
RTF (Microsoft Word互換) 'ファイル名.rtf'
Images: JPG and PNG 'ファイル名.png'
gt テーブルオブジェクト(R) 'gt

指定もかんたんで、例えばLatexコードを出力したいときは

msummary(regs, 'latex')

と打つだけでOK.
ファイルに保存したい場合は

msummary(regs, 'result.tex')  # Tex ファイル
msummary(regs, 'result.png') # PNGファイル

情報の編集

テーブル出力パッケージとして簡単にきれいな表を出力できることはわかったが、問題はどの程度柔軟に情報が編集できるかである。

見た目

まずテーブルの見た目については、modelsummaryはkableExtraやgtの柔軟性を応用して編集できるようになっているようだ。
私はgtはまだ使ったことがないが、テーブルを画像出力する時に力を発揮するらしい。
Rmarkdown/knitrで出力するケースではknitr::kableのエクステンションであるkableExtraが良い。
HTML用にもLatex用にも編集がしやすく、きれいなテーブルを出力できる。

htmlについてはgtパッケージが、latexについてはhtmlパッケージがデフォルトになっているらしいが、次のようなオプションで変更することも可能らしい。



例えばhtmlのテーブルをgtパッケージの関数で修飾してみると以下のようになる。

msummary(regs, output = 'html')%>%
  gt::tab_style(style = cell_fill(color = "lightcyan"), # 行に色を付ける
            locations = cells_body(rows = 3)) %>%
  gt::tab_style(style = cell_text(weight = "bold"),   #文字を太くする
            locations = cells_body(rows = 3))

f:id:keita43a:20200529203838p:plain

回帰分析表の要素

回帰分析表の要素については基本的な編集機能が実装されているようだ。

標準誤差・t統計量・p値・信頼区間

回帰分析表の係数の推定値の下に表示される統計量。経済学では標準誤差を示すのが一般的だが、他の統計量もstatisticオプションを変更することで示すことができる。

msummary(regs, statistic = 'std.error') # 標準誤差 (デフォルト)
msummary(regs, statistic = 'p.value')   # p値
msummary(regs, statistic = 'statistic') # t統計量
msummary(regs, statistic = 'conf.int', conf_level = .99) # 99% 信頼区間

f:id:keita43a:20200529203909p:plain
99%信頼区間を表示した場合

表のタイトルとメモ

表のタイトルもかんたん。titleオプションで指定します。

メモも簡単にnotesオプションで追加できるし、複数のメモをリストで表示できる。

msummary(regs,
         title = "表１：mtcarsの結果",          # タイトル
         notes = list('メモ1. 南無妙法蓮華経', # メモ
                      'メモ２. アーメン'))

f:id:keita43a:20200529203948p:plain

変数名の変更

変数名の変更は、coef_mapで行う。
予め変更前後の変数名の対応表をベクトルにしておくと便利そうだ。
新しく入るほうが右に来るらしい。(dplyrのrenameと逆なので注意）

# 変数名対応のベクトル
var_nam = c("mpg" = "Mile per gallon", "(Intercept)" = "Constant","cyl" = "Cylinder", "disp" = "Displacement",
            "hp" = "Horse Power")

msummary(regs, coef_map = var_nam)

f:id:keita43a:20200529195311p:plain

モデルフィットなどの統計量

msummary()ではデフォルトでR二乗などのモデルフィットの統計量が表示されている。
不必要な統計量はgof_omitオプションで非表示にできる。

msummary(regs, gof_omit = "R2|AIC|BIC")

f:id:keita43a:20200529200117p:plain

統計量のリストはgof_mapで見ることができる。

> modelsummary::gof_map
# A tibble: 14 x 4
   raw           clean         fmt   omit 
   <chr>         <chr>         <chr> <lgl>
 1 nobs          Num.Obs.      %.0f  FALSE
 2 r.squared     R2            %.3f  FALSE
 3 adj.r.squared Adj.R2        %.3f  FALSE
 4 AIC           AIC           %.1f  FALSE
 5 BIC           BIC           %.1f  FALSE
 6 logLik        Log.Lik.      %.3f  FALSE
 7 deviance      Deviance      %.2f  TRUE 
 8 df.residual   DF Resid      %.0f  TRUE 
 9 df.null       DF Null       %.0f  TRUE 
10 sigma         Sigma         %.3f  TRUE 
11 statistic     Statistics    %.3f  TRUE 
12 p.value       p             %.3f  TRUE 
13 df            DF            %.0f  TRUE 
14 null.deviance Deviance Null %.2f  TRUE
P値の星

stargazerと違って、msummaryではp値の星はデフォルトではつかない。（最近は星つけないジャーナルも多い）
starオプションをTRUEにすると、デフォルト（10%, 5%, 1%)の基準で星をつけてくれる。
もちろんカスタムも可能。"*"を変更することで好きな文字にすることもできる。

msummary(regs, stars = TRUE) # デフォルト, `* p < 0.1, ** p < 0.05, *** p < 0.01`
msummary(regs, stars = c("*" = .05, "**" = .01, "***" = .001)) # カスタム

f:id:keita43a:20200529200933p:plain
カスタムの結果

小数点以下の調整

小数点以下の調整はfmtオプションで行う。sprintfのフォーマットをとるので、
例えば小数点第２位までを表示したければ以下のように書く。

msummary(regs, fmt = '%.2f') # 小数点以下第2位まで


f:id:keita43a:20200529201440p:plain

huxregと違って、observationの数はきちんと整数になっているのがいいですね。

ちなみに、科学的表記法にする場合はfをeに変える。

> # 例
> sprintf("%.3f", 0.5342)
[1] "0.534"
> sprintf("%.3e", 0.5342)
[1] "5.342e-01
行の追加

例えば固定効果を入れたかどうか(Yes/No)などを書いた行を追加したい場合、stargazerではそれ専用のオプションが合ったが、
msummaryではhuxregのようにカスタムされた行を追加することで対応するらしい。
しかし、huxregより追加は容易である。

add_rowsオプションで追加する。表と同じ行数の文字ベクトルを作成して、listにして追加する。
(2020年９月19日追記) Ver. 0.4.1 より add_rows はデータフレームしか受け付けなくなったので、文字ベクトルではなくデータフレームで情報を与える必要があります。

##row1 <- c('固定効果1', 'No', 'No')  # Ver 0.4.1.以前
##row2 <- c('固定効果2', "No","No") # Ver 0.4.1.以前

# 0.4.1.よりデータフレームで与える
# tibble::tribbleで行ごとにデータを与えてtibbleを作れる。

rows <- tribble(~term, ~`Model A`, ~`Model B`,
               'FE 1', 'No', 'No',
               'FE 2', "No","No")  # Sep 20, added. Now the add_rows accepts only data.frame

attr(rows, 'position') <- c(15,16) # position attributeでどこに挿入するか指定する。15-16行目に追加

## msummary(regs, add_rows = list(row1, row2))# Ver 0.4.1.以前
msummary(regs, add_rows = rows)

f:id:keita43a:20200529202212p:plain

また、上の表な手作業ではなく、モデルをカスタマイズして固定効果を表示する方法も提案されていますが、別記事でまとめたいと思います。

で、estimatrは使えるのか？

大変使いやすそうだし、broom::tidyの出力に基づいているので柔軟そうだ。

ではestimatr::lm_robustは使えるのか？
結論から言うと、使えます!

regs_rob <- list()
regs_rob[["Model A rob"]] <- lm_robust(mpg ~ cyl+disp, 
                                       se_type = "stata",
                                       fixed_effects = ~ vs + am,mtcars)
regs_rob[["Model B rob"]] <- lm_robust(mpg ~ cyl+disp+hp, 
                                       se_type = "stata",
                                       fixed_effects = ~ vs + am,mtcars)

msummary(regs_rob)

f:id:keita43a:20200529202751p:plain


自動で標準誤差のタイプまで表示してくれていますね。

まとめ

  • {modelsummary}はシンプルで使いやすい回帰分析表出力パッケージ
  • broom::tidyの結果を受け付けるので、最近開発された統計モデルパッケージ等にも対応
  • gtやkableExtraを使ってテーブルが修飾できるし、様々なフォーマットへの出力に対応 (html, latex, RTF, text, 画像ファイル)
  • 多くのオプションでテーブルの内容もカスタマイズ可能

まだ実際に論文書くのには使ってはいませんが、よさそうなので使ってみます！
使ったコードはここにおいてます。

github.com

【R】スクレイピング時のループで404エラーが出た時に止まるのを避ける

R データ分析

ちょっと必要になって、ウェブサイト上のテーブルをまとめてダウンロードするという作業が必要になった。

だいたいPythonとかでやるのが主流なようだが、私はR使いなのでRでやる方法を探しているとrvestというパッケージがよいらしい。

rvestの基本的な使い方に関しては、Shohei Doiさんのノートが大変参考になった。
shohei-doi.github.io


さて、問題は同じようなウェブサイトページにURLをループで取得したい時に、時々ページ自体が存在せず、404 Error Not Foundのエラーが出ることがある。

これが困るのは、そのエラーが出た時点でループが止まってしまうことだ。とりあえず404があれば無視して次に行ってほしいので、どうすればいいかを調べてみたことをメモっておく。

tryCatchを使う。

以下のコードが実際に私が書いたコードの一部である。

ある記事に番号が振られており、1から300まで回して中にある表を取得しようとしているのだが、
不規則に404 Not Foundがでてくるので、単純にループを回すとそこで止まってしまう。

まず、delayedAssignでdo.nextという変数にnextの関数を割り振る。nextはloop内で出てくると、次のループに飛ばしてくれる命令である。

tryCatchは、１つ目の引数を評価して、エラーなど異常な条件が出ると、それに対して対応を決められる関数である。ここではrvestのread_htmlの関数を評価しているので、これにエラーがでると次のerror =で対応を決めている。
ここでは、もしURLが404 Not Foundでエラーが出た場合、"URL Not Found, skipping"というメッセージを表示して、force(do.next)により、次のループに飛ぶように命令している。

library(tidyverse)
library(rvest)

yr = 2017  #extracting year

list_tables = list()  # make list to store the result

for(i in 1:300){ # j-meldings up to 300 per year? maybe more in other years...
  
 url <- paste0("https://www.fiskeridir.no/Yrkesfiske/Regelverk-og-reguleringer/J-meldinger/Utgaatte-J-meldinger/J-",i,"-",yr)


 delayedAssign("do.next", {next})
 
 # use tryCatch to avoid stopping by error reading 404 page. 
 html <- tryCatch(
   read_html(url),
   error=function(e) {print("URL Not Found, skipping")
     force(do.next)})

 # check if <head> and <body> are properly read.
 
 # read html file from the url
 html <- url %>% read_html

 # print the number of meldings 
 print(paste0("melding ",i," is successfully read"))
 
# 2.1. check the title of melding
 
 # if the read html is not about the cod fishery, skip to next loop
 html_title <- html %>%
   html_node("h1") %>% 
   html_text() %>%
   str_detect(pattern = "Forskrift om regulering av fisket etter torsk, hyse og sei nord for 62°N")
 
 if(!html_title){
   next
 }

固定効果とランダム効果：統計学と計量経済学での定義

データ分析 計量経済学

ある空間的な変数を持つデータを扱っているときに、どういう形で空間的な構造を扱うか悩んでいた。
空間計量経済学はまだきちんと勉強したわけではなく、実証的な論文もそんなに読んでいないので扱いがいまいちわからない。

自分の学部には、Rの空間統計パッケージの開発などで超有名な地理経済の先生がいるので、*1エスプレッソマシンの近くでコーヒーブレークしていたところを捕まえて、質問してみた。*2

ざっとどんな分析がしたいか聞いてみたが「とりあえずマルチレベルモデルでいいんじゃない？」とのこと。計量経済学では聞き慣れない名前なので、どういう分析か聞いてみたところ、「計量経済学で言ういわゆるランダム効果(Random effect)ってやつだな。」と答えてくれた。「でも、統計学でのランダム効果と計量経済学では少し定義が異なるから、これを読むといいよ」と言って教えてくれたのが、Journal of Statistical SoftwareのRのplmパッケージの記事だった。

plmパッケージはパネルデータを扱うパッケージなのだが、基本は計量経済学におけるパネルデータ分析を念頭に設計されている。しかし、統計計算のソフトウェアを扱うジャーナルにおいて、計量経済学以外の統計分野の読者も多いため、査読の過程で計量経済学における用語の使い方と統計における使い方の整理をするセクションを書くように指示されたという経緯があったらしい。その過程で追加されたセクションが7.2のSome false friendsという項目だそうだ。

混合モデルとマルチレベルモデル

統計学や統計を用いる科学分野において、混合モデルやマルチレベルモデルはいくつかの複数の名前があるようだ。
例を上げてみると

混合モデル
マルチレベルモデル
階層線形モデル
混合効果モデル
ハイブリッドモデル

など。。。一般にはこれらのモデルは同じものを違う名前で呼んでいるとされるらしいが、
微妙に異なるコメントが以下のスレッドではなされていた。

stats.stackexchange.com


要約すると、マルチレベルモデルや階層線形モデルでは「マルチレベル」や「階層」という言葉から、入れ子になっている構造(Nested)が強調される。入れ子になっている構造とは、たとえば県の中に市町村があり、その中に個人が属しているようなデータで、県レベルの効果・市町村レベルの効果をそれぞれ推定するようなケースである。しかし、混合モデルにおいては必ずしも入れ子になっている必要はない。おそらくマルチレベルと呼びながら入れ子になっていない構造のデータを扱っている場合もあるのではないだろうか。


混合モデルと計量経済学のパネルデータ

用語の違いと対応

計量経済学の固定効果とランダム効果

以上の議論を受けて、冒頭ではマルチレベルモデルに出会った経緯を紹介したが、より広い定義を持つ混合モデル(mixed model)と呼ぶことにする。

計量経済学が統計学と少し離れて発展してきた背景には、社会科学を扱う学問として、非実験データ（観察データ）を扱ってきたため、観察データ特有の問題に対処する必要があったことが挙げられる。
特に問題となってきたのが内生性である。つまり、経済理論を実証する形で推定される計量経済学では経済理論の記述する変数間の関係が統計的に実証できるかどうかに関心が置かれるため、変数Xが変数Yに与える因果関係にバイアスがないかどうかがおおきな関心となる。そこでのパネルデータを用いた分析がクロス・セクションデータの分析より秀でているとされる点は、個人や時間の効果という直接観察されないが、興味のある変数Xと相関している場合にバイアスがかかるため(Omitted Variable Bias)、個人や時間について複数の観察点がある場合にはその効果が推定されることでバイアスが避けられる（可能性がある）ということである。

計量経済学におけるパネルデータ分析で推定される個人や時間の効果は、固定効果(Fixed Effects)とランダム効果(Random Effects)に分けられる。式で書くとどちらも

y_{it} = a + x_{it}'\beta + v_{it}

となり、 y_{it}が従属変数、 x_{it}が説明変数、 v_{it}が誤差項であるが、そのうち v_{it} = \mu_{i} + \varepsilon_{it} と分解され、 \mu_{i}が個人の効果で、 \varepsilon_{it}が平均ゼロの分布を持つ誤差項である。

この個人の効果が、個人特有の切片として推定される場合は固定効果、誤差項の一部としてパラメータを持つ分布として推定される場合はランダム効果と呼ぶ。
固定効果として推定される場合は、個人の効果は変数Xと相関を持つが、ランダム効果の場合は相関がないという仮定のもとで推定される。

計量経済学的にとくに関心があるのは固定効果の場合である。つまり、個人特有のなにか（直接観察されない）が変数Yに影響を与えるが、Xとも相関しているのでXの効果にバイアスが生じる場合、
固定効果を含めることでこのバイアスが避けられるというアイデアである。

統計学(混合モデル)の固定効果とランダム効果


混合モデルとは、固定効果とランダム効果の両方を含む回帰分析モデルを指す。

混乱のもとは、この固定効果とランダム効果(変量効果)という専門用語が、統計学と計量経済学の間で異なる定義で、しかし似通った文脈で使われていることである。

混合モデルにおける固定効果とは、パラメータが定数であるモデル（つまり普通の回帰分析）であり、ランダム効果とはパラメータが平均ゼロの同時正規分布に基づいてランダムに変化するようなモデルである。

すなわち、計量経済学におけるランダム効果とは、混合モデルのランダム効果のうち特殊なバージョンであり、パラメータのうち切片だけが分布に基づくランダム性を持つという形になる。

計量経済学における固定効果にあたる名前は混合モデルにはないが、あえていうならグループや個人ごとのダミー変数を用いて推定するダミー変数最小二乗法(DVLS)と言える。

逆に、混合モデルのランダム効果は、計量経済学では最近ではランダム係数モデル(Random Coefficient)と言われるモデルがある。

また、混合モデルの固定効果とは、上述のように普通の線形モデルである。


推定方法の違い

厳密には計量経済学と統計学の違い、というわけではないかもしれないが、特にランダム効果モデルの推定方法は異なるようだ。

計量経済学では、経済理論に基づいた変数間の関係の実証に興味がある。そのため、推定方法として一般的なのが(係数の推定自体は)分布を特定しないOLSである。
このOLSから派生したFealsible GLSを用いることで、一般的なOLSの過程を満たさない（均一分散ではない）ランダム効果モデルの推定を行う。

推定したいモデルが

y_{it} = a + x_{it}'\beta + v_{it}

であり、上と同様に誤差項が v_{it} = \mu_{i} + \varepsilon_{it} と分解される場合、ランダム効果モデルでは、まず外生仮定を以下のように仮定する。

E(\varepsilon_{it} | {\bf x_{i}} \mu_{i}) \forall t

そして、個人効果の条件付き期待値がゼロだと仮定する。
E(\mu_{i} | {\bf x_{i}}) = E(\mu_{i}) = 0

さらに、GLSの推定の仮定としてフルランク仮定を置く。

\text{rank } E(X_{i}'\Omega X_{i}) = K

K は変数の数であり、 \Omega = E(v_{i}v_{i} ') は無制約分散推定量である。

これにより、

\Omega = \left( \begin{array}{ccc} E(v_{i1}^{2}) & E(v_{i1}v_{i2}) & \ldots & E(v_{i1}v_{iT})\\ E(v_{i2}v_{i1}) & E(v_{i2}^2) & \ldots & E(v_{i1}v_{iT})\\ \vdots & & & \vdots \\\\ E(v_{i1}v_{iT}) & & \ldots & E(v_{iT}^2)\\ \end{array} \right)

仮定を適用すると、 E(v_{it}^2) = E(\mu_{i}^2) + 2E(\mu_{i}\varepsilon_{it}) + E(\varepsilon^2_{it}) = \sigma_{\mu}^2 + \sigma_{\varepsilon}^2と E(v_{it}v_{is} = E(\mu_{i}^2) = \sigma_{\mu}^2となる。
このうち \sigma_{\mu}^2が、ランダム効果の分布の分散の推定値となる。

この辺りの詳しい説明はパパ・ウールドリッジの10.4.1.節が詳しい。*3

Econometric Analysis of Cross Section and Panel Data, second edition (The MIT Press)

Econometric Analysis of Cross Section and Panel Data, second edition (The MIT Press)

Amazon



一方で、統計学ではこのようなランダム効果モデルは分布を特定化した状態で最尤推定を行うのが一般的なようだ。
理論上は、誤差項の分布が正規分布だと仮定して推定を行えば推定値はGLSと一致するはずである。

まとめ

以下に違いを表にまとめてみる。

昨今とくに、データ分析の隆盛に伴って分野を超えた議論が盛んになっているので、こういった違いが存在することを注意しながら議論・勉強していきたい。

項目 統計学(混合モデル) 計量経済学 (パネルデータモデル)
ランダム効果
  • 分布を持つランダム係数 (説明変数にかかる係数)
  • y_{it} = \alpha + x_{it}'\beta_{i} + \mu_{i} + \varepsilon_{it}
  • 分布を持つランダム切片
  • y_{it} = \alpha + x_{it}'\beta + \mu_{i} + \varepsilon_{it}
固定効果
  • 定数の係数
  • y_{it} = \alpha + x_{it}'\beta + \varepsilon_{it}
  • レベルごとの切片
  • y_{it} = \alpha + x_{it}'\beta + \mu_{i} + \varepsilon_{it}
推定方法 最尤法 Feasible GLS

*1:全然関係ないが、これだけパッケージ開発とかしてるけど、この先生のキーボードの打ち方が人差し指だけで打ついわゆるオヤジ打ちだったのをみて、やり方なんて関係なくて結果出した人が強いということを再確認した。

*2:言うまでもなくコロナ禍が起こる少し前の話である。ノルウェーでは３月中旬以降ロックダウンでオフィスへアクセスできていない。。

*3:Introductory Econometricsの方がベイビー・ウールドリッジ

【ノルウェー】歯医者さんの気になるお値段

海外生活 ノルウェー

ノルウェーで初めて歯医者に行った。

よくアメリカで歯医者に行くと治療費がめっちゃ高いから日本で行っておくように！と言われるが、アメリカ在住時に詰め物が外れてしまって歯医者に行ったときは、75ドルのところを大学の保険効かせて50ドルにしてもらった覚えがある。

ノルウェーはそんなもんじゃない。社会福祉が充実していて医療費が安いと言われるノルウェーであるが、なぜか歯医者だけは保険が基本的には効かず、高い。*1


しかし、数週間前から少し右下の奥歯が痛むというか違和感があったので、知り合いで歯科の博士課程にいらっしゃる方に相談したところ、よい歯医者さんを紹介してもらった。
本当が大学病院の学生がトレーニングでやる歯医者が安いらしいのだが、コロナウイルスの影響で閉まっているらしい。


ノルウェーの歯科医

その歯医者さんもロックダウン中だったのですぐには予約が取れず、ようやく今日行ってきた。
日本に留学経験のある中国人の先生で、すこし日本語を話してくれる。
オスロには日本人の歯医者さんがいるらしいが、ベルゲンにはいないので、言語に自信がなければこちらの先生がいいかもしれない。
（といっても治療中は全部英語だったが）

雰囲気は日本の歯医者と変わらないが、とにかく機材が多い印象。

問診票を出したあとは、すぐにレントゲン室に撮影に連れて行かれて、頭の周りをぐるっと廻る機械で撮影した後はすぐに診察室へ。
おそらく違和感は親知らずが圧をかけてるせいだろうということで、診察室の中でカメラみたいなものでレントゲンを追加で撮影した。
診察台の上にモニターがついていて、撮られた写真がすぐ見れる。こういう機材も今は日本でも普通なのかな。

とりあえず診断としては、緊急性はないが抜いたほうがいいと思う、という話だった。
そのお値段はなんと一本あたり3000クローネ程度だけど、知り合いの紹介ということで4本抜くなら10000クローネにしてくれるらしい。
日本円なら今だと約10万円か。これでも政府からの補助が出ていてこの値段らしい。
これでも以前別の歯科に妻を連れていったときの見積もりに比べると良心的だ。(20万円近かった覚えがある）

日本だとどれぐらいするんだろうか？とその場では見当がつかなかったので、とりあえず抜くかどうかはもう少し考えたいと伝える。
その後チェックしてもらったら、前歯(側切歯)の裏に左右それぞれ虫歯があるとのことだったので、そちらは処置してもらうことに。
時間もあったせいか、日本みたいに細切れで何回も通院する感じではなく、その場で全部処置してくれた。
おもったより深かったらしく、麻酔までしてもらうことに。
結局全部で1時間ぐらいかけて処置してもらった。

気になる価格は

気になるお値段だが、以下のような感じ。ちなみに本日時点でのJPY/NOK は10.48円である

ノルウェー語 日本語 価格(NOK)
Undersøkelse inkludert 2 bilder 検査(レントゲン二枚含む 1180
panoramarøntgen パノラマレントゲン 375
Hvit fylling 虫歯治療と詰め物 (左側切歯) 870
Hvit fylling 虫歯治療と詰め物 (右側切歯) 870
Apikal bilde 追加写真？ 0
Trygderefusjon 保険払い戻し -435
合計 2860

一応親知らずについては保険が少し効くらしい。私は住民番号を持っているので、保険が適用された。

それでも、仕方ないとはいえ3万円近い出費は痛い・・・



他に医療に関する記事も記録してます。

keita43a.hatenablog.com

keita43a.hatenablog.com

keita43a.hatenablog.com

*1:18歳未満の子供は無料で、19〜20歳までは75％がカバーされる。なお歯列矯正は無料の対象外。