餡子付゛録゛

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

ゲーム理論で考える、じゃんけんの拡張の数値演算

ゲーム理論で考える、じゃんけんの拡張: ニュースの社会科学的な裏側で使ったコードです。通常のゲーム理論の数値演算で扱われる計算手順に習っていないと言うか、フィーリングで描いたので、何か勘違いがあるかも知れないです。

##########################
### 引き分けになる確率 ###
##########################

p_even <- function(n){
  if(1>n){
    return(0)
  }
  1-(2^n-2)/3^(n-1)
}


##########################
### 確定的確率質量関数計算 ###
##########################

mkdist <- function(p){
  d2b <- function(num, digits=8){
    r <- numeric(digits)
    for(i in (digits-1):0){
      if(2^i <= num){
        num <- num - 2^i
        r[digits - i] <- 1
      }
    }
    return(r)
  }
  
  V <- matrix(0, 2^length(p), length(p)) # 組み合わせ
  pV <- numeric(2^length(p)) # 組み合わせの発生確率
  sV <- numeric(2^length(p)) # 戦略1の数
  for(i in 1:(2^length(p))){
    V[i, ] <- d2b(i-1, digits=length(p))
    pV[i] <- prod((0==V[i, ])*(1-p) + (1==V[i, ])*p)
    sV[i] <- sum(V[i, ])
  }
  tapply(pV, sV, sum) # 戦略1の数ごとの生起確率
}


######################
### 勝利確率の計算 ###
######################

calc_p_win <- function(p){
  p_win_sum <- 0

  # ぐー/ちょき/ぱーを出した他の参加者数ごとに計算していく
  for(num_of_nq in 0:length(p_others)){
  
    # ぐー/ちょき/ぱーで勝つ確率
    if(0==num_of_nq){
      # 自分だけぐー/ちょき/ぱーならば、確実に勝利
      p_win_nq <- 1
    }else{
      # あいこにならなければ、半分の確率で勝利
      p_win_nq <- (1 - p_even(num_of_nq + 1))/2
    }
  
    # きゅーで勝つ確率
    if(2>num_of_nq){
      # 一人がぐー/ちょき/ぱーならば、確実に敗北
      p_win_q <- 0
    } else {
      # あいこになれば、勝利
      p_win_q <- p_even(num_of_nq)
    }

    # 勝利確率
    p_win <- (1-p)*p_win_nq + p*p_win_q

    # ぐー/ちょき/ぱーで負ける確率
    if(0==num_of_nq){
      # 自分だけぐー/ちょき/ぱーならば、確実に負けない
      p_lose_nq <- 0
    }else{
      # あいこになれば敗北
      # あいこ以外でも半分の確率で敗北
      p_lose_nq <- p_even(num_of_nq + 1) + (1 - p_even(num_of_nq + 1))/2
    }
  
    # きゅーで負ける確率
    if(0 == num_of_nq){
      # 全員がきゅーならば、引き分けで負けない
      p_lose_q <- 0
    }else if(1 == num_of_nq){
      # 一人がぐー/ちょき/ぱーならば、確実に敗北
      p_lose_q <- 1
    }else{
      # あいこ以外ならば、敗北
      p_lose_q <- 1 - p_even(num_of_nq)
    }

    # 敗北確率
    p_lose <- (1-p)*p_lose_nq + p*p_lose_q

    # num_of_nqが生じる確率
    s <- as.character(length(p_others) - num_of_nq)
    p_num_of_nq <- dist[s][[1]]
    if(is.na(p_num_of_nq)){
      p_num_of_nq <- 0
    }

    # 期待値調整をして勝利確率を合計する(引き分けは分母から除く)
    p_win_sum <- p_win_sum + p_num_of_nq * p_win/(p_win+p_lose)
  
    if(debug){
      print(sprintf("num_of_nq_of_others: %d (%.3f) p_win_nq: %f p_win_q: %f p_win: %f p_lose_nq: %f p_lose_q: %f p_lose: %f", num_of_nq, p_num_of_nq, p_win_nq, p_win_q, p_win, p_lose_nq, p_lose_q, p_lose))
    }
  }
  p_win_sum
}


####################
### 均衡値の計算 ###
####################

# 拡張じゃんけん参加人数
n <- 5

# 初期値
p_all <- runif(n)

# n*10回ぐらい回せば収束するであろうと言う粗雑な方針
debug <- FALSE
for(c in 0:(n*10)){
  # 最適化前の状態
  print(sprintf("%.5f", p_all))

  # 最適化を行なうプレイヤー
  i <- (c %% n) + 1

  # 最適化を行なうi以外のプレイヤーのpを固定して、
  # キューを選択する人数の確率質量関数を作る
  p_others <- p_all[-i]
  dist <- mkdist(p_others)

  # 勝利確率の最大化を行なう
  r_optimize <- optimize(calc_p_win, c(0, 1), maximum=TRUE)
  p_all[i] <- r_optimize$maximum
}

print(sprintf("%.5f", p_all))
sprintf("最適反応の平均値: %.5f", mean(p_all))


#############################################################
# p_allから負ける人数の期待値を計算し、通常じゃんけんと比較 #
#############################################################
dist <- mkdist(p_all)
expected_num_of_looser <- 0
compared <- (1-p_even(length(p_all)))*length(p_all)/2
for(num_of_nq in 0:length(p_all)){

  s <- as.character(length(p_all) - num_of_nq)
  p_num_of_nq <- dist[s][[1]]

  t <- 0
  if(1<num_of_nq){
    if(1==num_of_nq){
      # 一人だけぐー/ちょき/ぱー
      t <- length(p_all) - 1
    } else {
      p_e <- p_even(num_of_nq)
      t <- p_e*num_of_nq + (1-p_e)*(num_of_nq/2 + length(p_all) - num_of_nq)
    }
    expected_num_of_looser <- expected_num_of_looser + p_num_of_nq*t
  }
}
sprintf("負ける人数の期待値 通常じゃんけん: %.5f 拡張じゃんけん: %.5f",compared, expected_num_of_looser)

多重共線性を出してみよう

観測数100ぐらいでも、誤差項の分散など次第では多重共線性が推定結果に影響を与える例を作ってみました。

多重共線性が無いケース

set.seed(2103)
n <- 100
x1 <- 0.3 * runif(n, min=0, max=10) + rnorm(n, sd=1)
x2 <- 0.7 * runif(n, min=0, max=10) + rnorm(n, sd=1)
cor(x1, x2) # 説明変数間は低い相関係数
e <- rnorm(n, sd=10)
y <- 1 + 2*x1 + 3*x2 + e
summary(lm(y ~ x1 + x2))

真のモデルにx2が無い場合に、x1とx2で推定してもそれらしい推定結果になる。

y <- 1 + 2*x1 + e
summary(lm(y ~ x1 + x2))

多重共線性があるケース

係数の大きさが変化して、有意性が無くなるケース。誤差項の分散は無いケースと同じ程度。

set.seed(2103)
n <- 100
x <- runif(n, min=0, max=10)
x1 <- 0.3 * x + rnorm(n, sd=0.5)
x2 <- 0.7 * x + rnorm(n, sd=0.5)
cor(x1, x2) # 説明変数間の相関係数が高い
e <- rnorm(n, sd=10)
y <- 1 + 2*x1 + 3*x2 + e
summary(lm(y ~ x1 + x2))

真のモデルにx2が無い場合に、x1とx2で推定するとx1の有意性も消える。

y <- 1 + 2*x1 + e
summary(lm(y ~ x1 + x2))

拡張Dickey-Fuller検定を横断面方向に拡張する

同一のデータ生成プロセスを踏んでいる複数の系列 y1, y2, y3, … を、それぞれ個別に単位根検定にかけると、定常であるのに標本数が不足して単位根の棄却に失敗しがちになる一方、c(y1, y2, y3, …) と単純に接続して単位根検定した場合、接続面に断続が生じるため単位根を過って棄却する事になります。

代表的な単位根検定のDickey-Fuller検定と、それに系列相関やトレンド項の影響も考慮したADF 検定では、差分データを前期の水準データで回帰して、前期の水準データの係数の標準誤差をτ分布で評価します。具体的には単位根が疑われる系列 Y_tを、\epsilonを誤差項として、

 Y_t = \beta Y_{t-1} + \epsilon

と書き、以下のように差分 \Delta Y_tを説明する式に変形し、

 Y_t - Y_{t-1} = \beta Y_{t-1} - Y_{t-1} + \epsilon
 \Delta Y_t  = (\beta - 1) Y_{t-1} + \epsilon = \gamma Y_{t-1} + \epsilon

 \gamma = 0であるかを検定し、棄却されれば単位根が無いとします。ここで異なる系列でも同一のデータ生成プロセスであれば、前期の水準データ Y_{t-1}の係数 \gammaは同一となります。つまり、系列ごとに差分データと前期の水準データの組み合わせを作成後にデータを統合すれば、複数の系列を同時に単位根検定する事ができます。

一般に時系列方向と横断面方向にデータの広がりのあるパネルデータの場合、横断面方向にデータの広がりがあるため、単一系列の場合と異なり見せかけの相関が生じにくく、株価などの不確実性の高い高頻度取引などではなく、観察対象が低頻度に理屈を立てて行動した結果でああるため、単位根の存在は気にしない事が多い*1わけですが、たまたま稀な需要があったのでurcaパッケージのur.dfを拡張してみました。

1. ソースコード

urcaを呼び出したあと、関数ur.dfを上書きします。ちょっと長いですがデータ統合する部分以外は、元と同じです。ですから、元と同じオプションが使えます。元のが使いたくなった場合は、環境を指定するか、追加定義したグローバル領域のur.dfをrm(ur.df)して消してください。

library(urca)
ur.df <- function (..., type = c("none", "drift", "trend"), lags = 1, selectlags = c("Fixed", "AIC", "BIC"))
{
  lag <- as.integer(lags)
  if (lag < 0)
    stop("\nLags must be set to an non negative integer value.\n")
  lags <- lags + 1

  args <- list(...)
  len <- length(args)
  if(len <= 0)
    stop("\ny must be set.\n")
  tmp <- as.list(NULL)
  tmp[["n"]] <- 0

  for(i in 1:len){
    y <- args[[i]]
    y <- y[!is.na(y)]

    if (ncol(as.matrix(y)) > 1)
      stop("\ny is not a vector or univariate time series.\n")
    if (any(is.na(y)))
      stop("\nNAs in y.\n")

    y <- as.vector(y)
    z <- diff(y)
    n <- length(z)
    x <- embed(z, lags)
    z.diff <- x[, 1]
    z.lag.1 <- y[lags:n]

    tmp[["y"]] <- append(tmp[["y"]], y)
    tmp[["z"]] <- append(tmp[["z"]], z)
    tmp[["n"]] <- tmp[["n"]] + n
    tmp[["x"]] <- rbind(tmp[["x"]], x)
    tmp[["z.diff"]] <- append(tmp[["z.diff"]], z.diff)
    tmp[["z.lag.1"]] <- append(tmp[["z.lag.1"]], z.lag.1)
    tmp[["tt"]] <- append(tmp[["tt"]], lags:n)
  }

  y <- tmp[["y"]]
  z <- tmp[["z"]]
  n <- tmp[["n"]]
  x <- tmp[["x"]]
  z.diff <- tmp[["z.diff"]]
  z.lag.1 <- tmp[["z.lag.1"]]
  tt <- tmp[["tt"]]

  (n-length(z.lag.1)+1):n

  selectlags <- match.arg(selectlags)
  type <- match.arg(type)

  CALL <- match.call()
  DNAME <- deparse(substitute(y))
  x.name <- deparse(substitute(y))

  if (lags > 1) {
    if (selectlags != "Fixed") {
      critRes <- rep(NA, lags)
      for (i in 2:(lags)) {
        z.diff.lag = x[, 2:i]
        if (type == "none")
         result <- lm(z.diff ~ z.lag.1 - 1 + z.diff.lag)
        if (type == "drift")
         result <- lm(z.diff ~ z.lag.1 + 1 + z.diff.lag)
        if (type == "trend")
         result <- lm(z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
        critRes[i] <- AIC(result, k = switch(selectlags,
         AIC = 2, BIC = log(length(z.diff))))
      }
      lags <- which.min(critRes)
    }
    z.diff.lag = x[, 2:lags]
    if (type == "none") {
      result <- lm(z.diff ~ z.lag.1 - 1 + z.diff.lag)
      tau <- coef(summary(result))[1, 3]
      teststat <- as.matrix(tau)
      colnames(teststat) <- "tau1"
    }
    if (type == "drift") {
      result <- lm(z.diff ~ z.lag.1 + 1 + z.diff.lag)
      tau <- coef(summary(result))[2, 3]
      phi1.reg <- lm(z.diff ~ -1 + z.diff.lag)
      phi1 <- anova(phi1.reg, result)$F[2]
      teststat <- as.matrix(t(c(tau, phi1)))
      colnames(teststat) <- c("tau2", "phi1")
    }
    if (type == "trend") {
      result <- lm(z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
      tau <- coef(summary(result))[2, 3]
      phi2.reg <- lm(z.diff ~ -1 + z.diff.lag)
      phi3.reg <- lm(z.diff ~ z.diff.lag)
      phi2 <- anova(phi2.reg, result)$F[2]
      phi3 <- anova(phi3.reg, result)$F[2]
      teststat <- as.matrix(t(c(tau, phi2, phi3)))
      colnames(teststat) <- c("tau3", "phi2", "phi3")
    }
  }
  else {
    if (type == "none") {
      result <- lm(z.diff ~ z.lag.1 - 1)
      tau <- coef(summary(result))[1, 3]
      teststat <- as.matrix(tau)
      colnames(teststat) <- "tau1"
    }
    if (type == "drift") {
      result <- lm(z.diff ~ z.lag.1 + 1)
      phi1.reg <- lm(z.diff ~ -1)
      phi1 <- anova(phi1.reg, result)$F[2]
      tau <- coef(summary(result))[2, 3]
      teststat <- as.matrix(t(c(tau, phi1)))
      colnames(teststat) <- c("tau2", "phi1")
    }
    if (type == "trend") {
      result <- lm(z.diff ~ z.lag.1 + 1 + tt)
      phi2.reg <- lm(z.diff ~ -1)
      phi3.reg <- lm(z.diff ~ 1)
      phi2 <- anova(phi2.reg, result)$F[2]
      phi3 <- anova(phi3.reg, result)$F[2]
      tau <- coef(summary(result))[2, 3]
      teststat <- as.matrix(t(c(tau, phi2, phi3)))
      colnames(teststat) <- c("tau3", "phi2", "phi3")
    }
  }
  rownames(teststat) <- "statistic"
  testreg <- summary(result)
  res <- residuals(testreg)
  if (n < 25)
    rowselec <- 1
  if (25 <= n & n < 50)
    rowselec <- 2
  if (50 <= n & n < 100)
    rowselec <- 3
  if (100 <= n & n < 250)
    rowselec <- 4
  if (250 <= n & n < 500)
    rowselec <- 5
  if (n >= 500)
    rowselec <- 6
  if (type == "none") {
    cval.tau1 <- rbind(c(-2.66, -1.95, -1.6), c(-2.62, -1.95,
      -1.61), c(-2.6, -1.95, -1.61), c(-2.58, -1.95, -1.62),
      c(-2.58, -1.95, -1.62), c(-2.58, -1.95, -1.62))
    cvals <- t(cval.tau1[rowselec, ])
    testnames <- "tau1"
  }
  if (type == "drift") {
    cval.tau2 <- rbind(c(-3.75, -3, -2.63), c(-3.58, -2.93,
      -2.6), c(-3.51, -2.89, -2.58), c(-3.46, -2.88, -2.57),
      c(-3.44, -2.87, -2.57), c(-3.43, -2.86, -2.57))
    cval.phi1 <- rbind(c(7.88, 5.18, 4.12), c(7.06, 4.86,
      3.94), c(6.7, 4.71, 3.86), c(6.52, 4.63, 3.81), c(6.47,
      4.61, 3.79), c(6.43, 4.59, 3.78))
    cvals <- rbind(cval.tau2[rowselec, ], cval.phi1[rowselec,
      ])
    testnames <- c("tau2", "phi1")
  }
  if (type == "trend") {
    cval.tau3 <- rbind(c(-4.38, -3.6, -3.24), c(-4.15, -3.5,
      -3.18), c(-4.04, -3.45, -3.15), c(-3.99, -3.43, -3.13),
      c(-3.98, -3.42, -3.13), c(-3.96, -3.41, -3.12))
    cval.phi2 <- rbind(c(8.21, 5.68, 4.67), c(7.02, 5.13,
      4.31), c(6.5, 4.88, 4.16), c(6.22, 4.75, 4.07), c(6.15,
      4.71, 4.05), c(6.09, 4.68, 4.03))
    cval.phi3 <- rbind(c(10.61, 7.24, 5.91), c(9.31, 6.73,
      5.61), c(8.73, 6.49, 5.47), c(8.43, 6.49, 5.47),
      c(8.34, 6.3, 5.36), c(8.27, 6.25, 5.34))
    cvals <- rbind(cval.tau3[rowselec, ], cval.phi2[rowselec,
      ], cval.phi3[rowselec, ])
    testnames <- c("tau3", "phi2", "phi3")
  }
  colnames(cvals) <- c("1pct", "5pct", "10pct")
  rownames(cvals) <- testnames
  new("ur.df", y = y, model = type, cval = cvals, lags = lag,
    teststat = teststat, testreg = testreg, res = res, test.name = "Augmented Dickey-Fuller Test")
}

2. 使い方

多少はと言った程度ですが、個別に単位根検定をかけるよりは、系列をまとめて検定する方が正しくなります。

genS <- function(n=200, beta=1.0){
  x <- numeric(n)
  x[1] <- 0
  for(i in 2:n){
    x[i] <- beta*x[i-1] + rnorm(1)
  }
  x
}

set.seed(121) # 数字を固定しないと結果はマチマチ
y1 <- genS(30) # 非定常過程を生成
y2 <- genS(30)
y3 <- genS(30)
summary(ur.df(c(y1, y2, y3))) # 単位根を過って棄却する
summary(ur.df(y1, y2, y3)) # 単位根を棄却しない

set.seed(100) # 数字を固定しないと結果はマチマチ
y4 <- genS(30, 0.9) # 定常過程を生成
y5 <- genS(30, 0.9)
y6 <- genS(30, 0.9)
summary(ur.df(y4)) # 単位根を棄却しない
summary(ur.df(y5)) # 単位根を棄却しない
summary(ur.df(y6)) # 単位根を棄却しない
summary(ur.df(y4, y5, y6)) # 単位根を棄却する

*1: 同時性、不均一分散、系列相関はとても気にされる。

打ち切りリセット付単位根過程のAR(1)での推定

ちょっと特異なデータを、そこそこ正しく推定するための検討用です。結論は、AR(1)にしておけば無問題でした。

### データ生成 ###
genSample <- function(n=250){
  ex <- rnorm(n)
  x <- rep(NA, length = n)
  x[1] <- ex[1] * 0.2
  for (j in 2:n) {
# 係数0.98で微妙に単位根ではないですが、250だとほぼ単位根過程と言うことで
    x[j] <- x[(j - 1)] * 0.98 + ex[j] * 0.2
  }
  x <- x + 5
  
  ey<-rnorm(n)
  y<-rep(NA,length=n)
  y[1]<-60+ey[1]*5
  
  d <- numeric(n)
  dc <- 0
  
  g <- numeric(n)
  gc <- 0
  
  chg <- numeric(n)
  
  for(i in 2:n){
    if(y[(i-1)]>20){
      y[i]<-y[(i-1)]-2+ey[i]*5
      dc <- dc + 1
      chg[i] <- 0
    }
    else{
      y[i]<-60+ey[i]*5
      dc <- 0
      gc <- gc + 1
      chg[i] <- 1
    }
    d[i] <- dc
    g[i] <- gc
  }
  data.frame(
    y2 = - y[-1],
    y1 = y[-length(y)],
    x = x[-1],
    d = d[-1],
    g = g[-1],
    chg = chg[-1]
  )
}

### m回施行して傾向を見る ###
n <- 250
m <- 3000
Pr_LVL <- numeric(m)
Pr_AR1 <- numeric(m)
for(c in 1:m){
  df1 <- genSample(n)

  r1 <- lm(y2 ~ x + factor(g) + d, data=df1)
  Pr_LVL[c] <- coef(summary(r1))[,4][2]

  r2 <- lm(y2 ~ y1 + x + factor(g) + chg, data=df1)
  Pr_AR1[c] <- coef(summary(r2))[,4][3]
}

hist(Pr_LVL, breaks=10)
hist(Pr_AR1, breaks=10)

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

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

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)
e1 <- exp(m %*% matrix(coef(r_ols), l, 1))
e2 <- exp(sum(coef(r_ols) * c(1, 0, 0, 0, 0)))
PCT_0to1 <- e1/(1+e1) - e2/(1+e2)
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){
  # 高低2種類の労働生産性
  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))