餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

オッズを使ってロジスティック回帰してみる

ロジット・モデルと言うと個票に応用する二項選択モデルを最尤法などで推定することを思い浮かべる人が多いと思いますが、都道府県ごとの率など集計データに一般線形回帰をかける事で推定する事もできます。観測数不明の集計データしか取れない場合も多いので、こちらの方が重宝するかも知れません。試しに使ってみましょう。

1. ロジット・モデルの線形化

ロジットモデルは以下のような累積分布関数から事象の発生確率pを定めます。

 p = \frac{ e^{x \beta} }{ 1 + e^{x \beta} }

xは条件を表す説明変数のベクトル、\betaはその係数です。 x \betaの値に関わらずpは0から1の間に収まるので、1超や負のあり得ない確率は出てこなくなるのが利点です。
このときオッズ p/(1-p)を整理すると、以下のようになります。

 \frac{p}{1-p} = e^{x \beta}

対数を取ると以下のように、線形モデルになります。

 \log \frac{p}{1-p} = x \beta = \alpha + \beta_1 x_1 + \beta_2 x_2 + \cdots

これに誤差項がつくと仮定して、重回帰をかけます。

2. 試しに推定してみる

手ごろなデータセットが『ブール代数分析による社会的カテゴリーの研究:「日本人」カテゴリー認識の分析*1の中にあったので拝借してきて、推定してみます。日本国籍を持ち社会学関係科目を履修している関西の大学生205名が、日本国籍・日本在住・日本人血統・日本語話者の4属性から日本人か否かをどう判別するかアンケートを取ったもので、サンプリング・バイアスと血統属性のトートロジー感(親が日本国籍?)が気になるデータですが、推定自体は問題なくできます。

### データセット ###
# N:日本国籍,R:日本在住,B:血統,L:日本語話者,ANS:日本人とする率
df1 <- data.frame(
  N = c(1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0),
  R = c(1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,0),
  B = c(1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0),
  L = c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0),
  ANS = c(0.995,0.922,0.922,0.659,0.995,0.756,0.82,0.502,0.639,0.4,0.278,0.059,0.493,0.18,0.049,0.015))

# オッズを計算する
attach(df1)
df1$ODDS <- df1$ANS/(1-df1$ANS)
detach(df1)

# 回帰分析をかける
r_ols <- lm(log(ODDS) ~ N + R + B + L,data=df1)
summary(r_ols) # 結果確認

指数をとらないと確率にならないので、推定結果だけでは分かりづらいですね。
親は日本人では無いものの、日本国籍がある日本語話者が日本人認定される確率を計算すると以下になります。

e <- exp(sum(coef(r_ols) * c(1, 1, 0, 0, 1)))
e/(1+e)

全ての説明変数が0を基準に限界効果を計算してみましょう。

e <- exp(sum(coef(r_ols) * c(1, 0, 0, 0, 0))) # c(...)の中の最初の1は切片項のため
coef(r_ols)[2:length(coef(r_ols))] * e/(1+e)^2 # e/(1+e)^2が密度関数になる

説明変数が不連続なので非直観的な数字になります。
これでは解釈しづらいので、0から1に変わったときの確率の計算もします。

l <- length(coef(r_ols))
m <- matrix(1, l-1, l)
m[,2:l] <- diag(l-1)
e <- exp(m %*% matrix(coef(r_ols), l, 1))
PCT_0to1 <- e/(1+e) - exp(sum(coef(r_ols) * c(1, 0, 0, 0, 0)))
rownames(PCT_0to1) <- c("N", "R", "B", "L")
PCT_0to1

3. 多値選択モデル(Multinomial Logit Model)への対応

3値以上の多値選択モデル(Multinomial Logit Model)でも、複数の重回帰を行なう事で推定値が出ます。
選択肢の数を J、選ばれた選択肢を mとして累積密度関数は、

 p_m = \frac{e^{x \beta_m}}{1 + \sum^{J-1}_{j=1} e^{x \beta_j} }

 p_0 = 1 - \sum^{J-1}_{j=1} p_j = \frac{1}{1 + \sum^{J-1}_{j=1} e^{x \beta_j} }

となります。ここで p_0は、直接推定されない参照カテゴリーが選ばれる確率です。 \betaが選択肢ごとに変わることに注意してください。
オッズは、

 \frac{p_m}{p_0} = e^{x \beta_m}

です。この式の対数をとって、 m=1,\cdots,Jの分だけ重回帰をかけると推定が終わりです。

となります。

4. 欠点

集計カテゴリーで差が出ないデータだと、集計することで説明変数と被説明変数の差異が小さくなります。この場合、検出力が露骨に落ちます。個票が使える場合は、なるべく個票を使いましょう。

おまけ:力技的に元データを復元

実はこのデータは、全ての説明変数の組み合わせの回答者数がそれぞれ205だと分かっているので、観測数16の集計データから、観測数16×206=3280の質的データに変換した上でglm関数で推定する事も可能です。

n <- 205
m <- 2^4 # 説明変数の数が4
df2 <- data.frame(
  N = numeric(n*m),
  R = numeric(n*m),
  B = numeric(n*m),
  L = numeric(n*m),
  ANS = numeric(n*m)
)
attach(df1)
for(i in 1:m){
  n1 <- round(n*ANS[i], 0)
  n2 <- n-n1
  b <- (i-1)*n + 1
  e <- i*n
  df2$N[b:e] <- rep(N[i], n)
  df2$R[b:e] <- rep(R[i], n)
  df2$B[b:e] <- rep(B[i], n)
  df2$L[b:e] <- rep(L[i], n)
  df2$ANS[b:e] <- c(rep(1, n1), rep(0, n2))
}
detach(df1)
r_glm <- glm(ANS ~ N + R + B + L, family = binomial(link = "logit"), data=df2)
summary(r_glm)

この推定結果も被説明変数をオッズとした推定結果になっています。結果はだいたい一緒ですが、微妙に違います。

*1: この論文ではブール代数分析と言って論理式と言うかベン図のようなもので日本人判定をしている事を仮定に分析しています。

観察された相対リスクが全て交絡バイアスによるものだとしても、未測定交絡因子の相対リスクはそんなに高くなくても良いかも知れない

受動喫煙防止法について論点整理①:受動喫煙による健康リスク・死亡者数の推定はどのくらい信用できるか?」と言うブログのエントリーで、受動喫煙の相対リスク1.22倍が未測定交絡因子によって引き起こされたバイアスだとして、未測定交絡因子の相対リスクは何倍ぐらい無いといけないかと言う議論をしていて、ノンパラメトリックな分析方法を使って最低でも1.74としているのですが、参照している分析手法の論文の記述をよく見ると、もっと低い値でも済みそうなので*1数値例を作ってみました。

set.seed(5301929)

### データ生成 ###

# 観測数
n <- 300

# 観察不能な1か0の未測定交絡因子
X <- c(rep(0, n/2), rep(1, n/2))

# 90%の確率で未測定交絡因子と一致する観察可能因子
CV <- 0.9
E <- c(ifelse(runif(n/2) < CV, 0, 1), ifelse(runif(n/2) < CV, 1, 0))

# ロジット関数
plogit <- function(p){
  exp(p)/(1 + exp(p))
}

# ロジットモデルにおける定数項と未測定交絡因子Xの係数
CONS <- -1
COEF_X <- 0.5

# イベント発生の有無を示す1か0かの被説明変数
D <- 0 + (runif(n) < plogit(CONS + COEF_X*X))*1

sprintf("未測定交絡因子Xの相対リスク: %.4f", plogit(CONS + COEF_X)/plogit(CONS + 0))
sprintf("未測定交絡因子Xと観察可能因子Eの相関係数: %.4f", cor(X, E))


# 観察可能因子と未測定交絡因子のリスク比RR_EXを計算する(注意:論文ではRR_{EU})
r_xtabs <- xtabs(~X + E, data=data.frame(D = D, X = X, E = E))
RR_EX <- max(r_xtabs[4]/r_xtabs[3], r_xtabs[2]/r_xtabs[1]) # = max(P(X=1|E=1)/P(X=1|E=0), P(X=0|E=1)/P(X=0|E=0))
sprintf("観察可能因子Eと未測定交絡因子Xのリスク比: %.4f", RR_EX)


### 説明変数をEで、ロジット関数の推定をする ###
r_E <- glm(D ~ E, family = binomial("logit"))
sprintf("観察可能因子Eの観察された相対リスク: %.4f", plogit(coef(r_E)[1] + coef(r_E)[2])/plogit(coef(r_E)[1]))

(追記:nを30000、CONSをlog(0.0038/(1-0.0038))、COEF_Xを0.32ぐらいにすると、実際の分析データに近い設定になります。)

未測定交絡因子の相対リスクは1.40倍で、観察可能因子Eの観察された相対リスクは1.35倍になりました。未測定交絡因子と観察可能因子の強い相関を仮定したら必要となる相対リスクは低いものの、そうでないと高くなると言う話のようです。

*1: the observed relative risk of cigarette smoking on lung cancer was RR^{obs}_{ED}=10.73 ... The work of Cornfield et al. showed that for a binary unmeasured confounder to completely explain away the observed relative risk, both the exposure-confounder relative risk and the confounder-outcome relative risk would have to be at least 10.73

労働市場シグナリング仮説と外部効果の推定値のシミュレーション

労働市場シグナリング仮説が成立していても、計量分析で推定される外部効果は負になるとは限らない事を確認するための、簡単なシミュレーションです。

#
# 労働市場シグナリング仮説に基づく賃金計算関数
# edu_h_s: 大卒比率
# return: 大卒賃金, 高卒賃金
#
wage <- function(edu_h_s){
  # 大卒と高卒の労働生産性
  lp_h <- 1.0
  lp_l <- 0.5

  # 労働生産性の高低の比率
  lp_h_s <- 0.5
  lp_l_s <- 0.5

  # 賃金=平均労働生産性
  if(lp_h_s > edu_h_s){
    # 高労働生産性の人の一部しか進学しないケース
    w_h <- lp_h
    w_l <- (lp_h*(lp_h_s-edu_h_s) + lp_l*lp_l_s)/(lp_h_s-edu_h_s+lp_l_s)
  }else{
    # 高労働生産性の人の全員と低労働生産性の人の一部が進学するケース
    w_h <- (lp_h*lp_h_s + lp_l*(edu_h_s-lp_h_s))/edu_h_s
    w_l <- lp_l
  }

  return(c(w_h, w_l))
}

set.seed(201705)

# 50地区10期のデータを生成する
noa <- 50 # Number of Areas
not <- 10 # Number of Spans

n <- noa*not
df <- data.frame(t=numeric(n), edu_h_s=numeric(n), wage_lp=numeric(n))
for(t in 1:not){
  b <- 1 + noa*(t-1)
  e <- b + noa - 1
  # 時点を定める
  df$t[b:e] <- t

  # まず、毎年、増加する、全体の大学進学率を定める
  edu_h_a <- 0.15 + 0.3*t/not
  # 各地域の大卒人口を確率的に設定
  edu_h_pop <- rnorm(noa - 1, mean=edu_h_a, sd=0.01)
  # 最後の地域を調整し、全体の大学進学率にあわせる
  edu_h_pop <- c(edu_h_pop, noa*edu_h_a - sum(edu_h_pop))

  # 全体の高卒比率は1-大学進学率
  edu_l_a <- 1 - edu_h_a
  # 各地域の高卒人口を確率的に設定
  edu_l_pop <- rnorm(noa - 1, mean=edu_l_a, sd=0.01)
  # 最後の地域を調整し、全体の高卒比率にあわせる
  edu_l_pop <- c(edu_l_pop, noa*edu_l_a - sum(edu_l_pop))

  # 各地域の大卒比率を計算する
  df$edu_h_s[b:e] <- edu_h_pop/(edu_h_pop + edu_l_pop)

  # 全体の教育水準の平均値から賃金を計算
  # 大卒と高卒が移動しても、それぞれの地域の大卒職場と高卒職場の平均労働生産性の期待値には変化が無いので、全体で計算している。生産性の高い高卒の移動が霍乱要因になるが、無数の人間が移動することで期待値に等しくなっていると解釈。
  df$wage_lp[b:e] <- wage(mean(df$edu_h_s[b:e]))[2]
}

#
# 回帰分析
# ・ミンサー・アプローチに習って対数化。しなくても同様の結果
# ・被説明変数は高卒賃金なので、大卒ダミーはつけていない
#

# 時系列ダミーが無いと、外部不経済が観測される
summary(lm(wage_lp ~ edu_h_s, data=df))

# 時系列ダミーがあると、外部不経済が観測されない
summary(lm(wage_lp ~ edu_h_s + t, data=df))

最適化レベルで結果が変わってしまうC言語のダメなコード

前世紀に書いたgccで-O1と-O3で実行結果が異なるコードです。インライン展開が作用しました。今のバージョンのgccでどうなるかは確認していません。なお、警告は出ると思いますし、処理系によってはコンパイルが通らない事もあると思います。

#include<stdio.h>

int *ptr;

void test(void){
  int a;
  a = 10;
  ptr = &a;
}

void test2(void){
  int b;
  b = 23;
}

main(){
  test();
  printf("Result1: %d\n",*ptr);
  test();
  test2();
  printf("Result2: %d\n",*ptr);
}

Rの関数をラッピング

Rは対話型インターフェイスなので、長々と関数に引数を指定するのは手間隙なので、グローバル・オプションを設定したくなります。しかし、発展途上のライブラリが多いため、グローバル・オプションをつけると上手く動作しなくなったり、そもそもグローバル・オプションを参照しないライブラリがあったりします。こういう都合の悪いときに、やや強引にデフォルト値を設定する細工を考えてみました。

R言語では既に宣言されている関数や変数を再定義する事ができ、さらに再定義した関数内からもともと宣言されている関数を呼び出す事ができます。これを使って、read.tableのfileEncodingにデフォルト値を設定してみましょう。

read.table <- function(...){
  print("ラッパーを通しています(理解したらここは消してください)")
  # 可変長引数をリストに展開
  args <- list(...)
  # fileEncodingが存在しなければ、UTF-8をセット
  if(!("fileEncoding" %in% names(args))){
    args <- append(args, list("fileEncoding"="UTF-8"))
  }
  # ホンモノのread.tableを呼び出す
  # package:utils はホンモノがあるEnvironment
  do.call(get("read.table", env=as.environment("package:utils")), args)
}

これを .Rprofile *1に書いておけば、options(encoding="UTF-8")なしにデフォルトでUTF-8指定ができます。write.tableでも、他の関数でも同様にできます。なお、ラッピングする関数の環境(Environment)が分からない場合は、find("関数名")としてください。
しかし、こういうコーディング、MS-DOS時代の割り込み処理を思い出しますね!

*1: Windowsの場合は、%USERPROFILE%\Documents\.Rprofile になります。

Rの文字コードのデフォルト値の設定

Windows版のRはテキストファイルの入出力にCP932と言うかShiftJISを使うわけですが、LinuxMacintosh OS XWindowsを併用して使っている人は、UTF-8をデフォルトにしたいかも知れません。read.tableのfileEncodingオプションをいちいち指定するのは面倒です。
解決方法は単純で、起動時に読み込む .Rprofile *1に、options(encoding="UTF-8") を書いておけば、UTF-8をデフォルトにします。設定値を確認したい場合は options()$encoding を表示してください。
なお、オブジェクトをそのまま保存するsave()とload()はコード変換が行なわれないので、使わない方が無難です。読み書き速度が問題になる場合は、readr パッケージを使うと改善できます。

*1: Windowsの場合は、%USERPROFILE%\Documents\.Rprofile になります

マクロ経済学の動的計画法説明用コードのMatlabからRへの移植

f:id:uncorrelated:20161224234007p:plain

一橋大学の阿部教授が院生向けのレクチャー・ノートをアップロードしていたのを見つけ、拝読して勉強させて頂いていたのですが、説明用のコード(P.16--20)がMatlabだったのでRに移植してみました。プロプライエタリMatlab嫌いな人に役立つかも知れないので、公開します。

単に移植するだけではなく、ループ周りを中心に無駄を省いています。理由は、行列にしておいた方が変数の中身を説明しやすいのと、コード全体の見通しが良くなるのと、速度面で有利になるからです。資本量の区切り幅inK=0.4ぐらいにして、行列の中身を表示させると、何を演算しているのか理解しやすいと思います。

alpha <- 0.25 # production parameter
beta <- 0.9 # subjective discount factor
gamma <- -2 # preference parameter
delta <- 1 # 1 -depreciation rate
minK <- 0.2 # minimum value of the capital grid
maxK <- 1.8 # maximum value of the capital grid
inK <- 0.001 # size of capital grid increments
nK <- ceiling((maxK-minK)/inK + 1) # number of grid points

# 横(列)の当期資本と、縦(行)の来期資本で定まる効用を表す行列
U <- matrix(rep(-Inf, nK^2), nK, nK)

K <- t(replicate(nK, ((1:nK)-1)*inK + minK)) # 当期資本。nK×nKに拡張。当期資本は列ごとに変えているのに注意。
K_prime <- replicate(nK, ((1:nK)-1)*inK + minK) # 来期資本k_{t+1}。nK×nKに拡張。来期資本は行ごとに変えているのに注意。
I <- K_prime - delta*K # 投資量。横(列)の当期資本と、縦(行)の来期資本で定まる。
C <- ((1-beta)/(alpha*beta))*(K^(alpha)) - I # (79)式を(80)式に代入した式から消費量を計算
U[C>0] <- ((C^(gamma+1))/(gamma+1))[C>0] # 効用を計算して行列に格納。消費がゼロ以下の場合は、初期値-Infを保存している

#
# whileループのための変数を初期化
#
V <- rep(0, nK) # value functionの代わりになるベクトル
Decis <- rep(0, nK) # Vに対応する来期資本を表すグリッド位置
iter <- 0 # ループ回数
metric <- 10 # 更新前後のvalue functionの差分の最大要素の値が入る。初期値は終了条件より大きければ、何でも良い。

while(metric > 0.001){
  # (82)式の最大化問題に対応する来期資本を探す
  # beta*Vは1列分のベクターだが、自動的にnK列分に複製拡大して足されるのに注意
  r <- U + beta*V

  # それぞれの列の最大値を選んで、value functionとする。
  # 当期資本に対応する列の中の最大値は、最大化問題の解となる来期資本
  tmpV <- apply(r, 2, max) # value functionの値
  tmpDecis <- apply(r, 2, which.max) # value functionに対応する来期資本を表すグリッド位置

  metric <- max(abs(V - tmpV)) # 終了条件と比較する更新によって生じた差分を計算

  V <- tmpV;
  Decis <- tmpDecis;
  iter <- iter + 1 # ループ回数をカウント

# 計算状況を表示
  vfor8 <- V[1] # value of V for k=kmin=0.2
  ufor8 <- (Decis[1] - 1)*inK + minK # value of control for k=0.2
  print(sprintf("ループ回数:%d, VF変化量:%.4f, vfor8:%.4f, ufor8:%.4f", iter, metric, vfor8, ufor8))
}

# policy functionと対応する各期の効用を計算
policy <- (Decis - 1)*inK + minK # policy function
U_op <- rep(-Inf, nK) # Uity under the optimal policy
K <- ((1:nK)-1)*inK + minK
I <- policy - delta*K;
C <- (1-beta)/(alpha*beta)*(K^alpha) - I
U_op[C>0] <- ((C^(gamma+1))/(gamma+1))[C>0]

# 各期の効用の割引価値を計算
betam <- beta*rep(1, nK)
value <- U_op/(rep(1, nK)-betam)

# 結果をプロット
x_at <- seq(1, length(value), length.out=17)
x_labels <- sprintf("%.2f", seq(minK, maxK, length.out=length(x_at)))
y_unit <- 5
ylim <- c(floor(min(value)/y_unit)*y_unit, ceiling(max(value)/y_unit)*y_unit)
y_at <- sprintf("%.2f", seq(ylim[1], ylim[2], length.out=5))
plot(value, main="DETERMINISTIC GROWTH MODEL: VALUE FUNCTION", xlab="Capital", ylab="Present Value", axes=FALSE, ylim=ylim)
axis(1, at=x_at, labels=x_labels, las=2)
axis(2, at=y_at)

線形近似の方は省略しました。