餡子付゛録゛

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

使っている名詞で2つの文の相関係数を取ってみる

Twitterで、ぬ氏が「二つの日本語の単語列(10-20文字ぐらい)があって、なるべく共通の単語が多いとか意味合い的に近いものを自動でマッチしたいときって何かいい方法ありますか?」と言うお題を出していたので、お気楽な方法を考えてみました。

わかち書きフレームワークのMeCabを利用して、2つの文を分解、単語1ダミー、単語2ダミー・・・と言う風に文Aと文Bのパラメーターをそれぞれつくって、相関係数を出して見ます。共通の単語が多ければ高い相関係数になるので、それで何かマッチング・アルゴリズムにかければ目的を達することができるでしょう。

説明すると長いわけですが、MeCabをインストールしてパスを通してコマンドラインで動くようにしてから、RMeCabをインストールして以下のように走らせたら終わりです。簡単ですね。

library(RMeCab)

# 比較する文を用意する
text <- c("牧瀬紅莉栖はオレの嫁", "牧瀬紅莉栖はツンデレ", "紅莉栖が生きている")

getWords <- function(text){
  res <- RMeCabC(text)
  value <- sapply(res, function(l){
    l[[1]]
  })
  type <- sapply(res, function(l){
    names(l[1])
  })
  value[type=="名詞"]
}

makeDummies <- function(words_all, words){
  0 + (words_all %in% words)*1
}


# 文ごとに名詞リストをつくる
words <- list()
for(i in 1:length(text)){
  words[[i]] <- getWords(text[i])
}

# 重複する単語は消して全体リストを作る
words_all <- unique(sort(unlist(words)))

# 文ごとに単語ダミーをつくる
words_dummies <- list()
for(i in 1:length(text)){
  words_dummies[[i]] <- makeDummies(words_all, words[[i]])
}

# 行列にまとめる
m <- matrix(unlist(words_dummies), length(words_all))
colnames(m) <- sprintf("text%d", 1:length(text))
rownames(m) <- words_all

# 相関係数を計算
cor(m)

ユーザー辞書を追加した、手元の環境だとこんな感じになりました。なお、ほとんどフィーリングで書いたので、RMeCabのコマンドを確認したら無駄が省ける可能性は大です。
f:id:uncorrelated:20200528183041p:plain

それでもすぐに、クラスター分析に持っていくことができます。

colnames(m) <- text
r_cluster <- hclust(dist(cor(m)), method="ward.D")
plot(r_cluster, main="MeCabで短文をクラスター分析", sub="", xlab="", ylab="", axes=FALSE)

f:id:uncorrelated:20200529112003p:plain

もちろん、具体的に何をするかは見えないので、役に立つかは謎です。文字列が長い場合は、よくあるテキストマイニング手法を真似しましょう。

同義語を揃える

使い道次第ですが、文書中に同義語(e.g. 社員募集=クルー募集)が多数考えられる場合は、同義語テーブルをつくって、一つの単語に揃えるほうが適切になるでしょう。

# 同義語テーブル(ベクターで用意/多い場合はファイルから読み込み推奨)
sv <- c("紅莉栖", "クリスティーナ", "助手", "セレブセブンティーン", "セレセブ", "蘇りし者", "ザ・ゾンビ")

# 最初の単語に揃える関数を用意(上にあわせ場合、MeCabに辞書登録している単語のみ有効)
l <- list()
for(i in 2:length(sv)){
  l[sv[i]] <- sv[1]
}

synonym <- function(words) sapply(words, function(word){
  ifelse(is.null(l[[word]]), word, l[[word]])
})

synonym(c("紅莉栖", "クリスティーナ", "セレセブ", "岡部"))

こんな感じで置換されます。

f:id:uncorrelated:20200529111723p:plain

MeCabのユーザー辞書

ユーザー辞書作成につかったcsvファイルの中身は以下です。

紅莉栖,1291,1291,5000,名詞,固有名詞,人名,名,*,*,くりす,クリス,クリス
クリスティーナ,1291,1291,5000,名詞,固有名詞,人名,名,*,*,くりすてぃーな,クリスティーナ,クリスティー
ツンデレ,1285,1285,1000,名詞,一般,*,*,*,*,*,ツンデレ,ツンデレ
セレブセブンティーン,1285,1285,5000,名詞,一般,*,*,*,*,*,セレブセブンティーン,セレブセブンティー
セレセブ,1285,1285,5000,名詞,一般,*,*,*,*,*,セレセブ,セレセブ
蘇りし者,1285,1285,5000,名詞,一般,*,*,*,*,*,ヨミガエリシモノ,ヨミガエリシモノ
ザ・ゾンビ,1285,1285,5000,名詞,一般,*,*,*,*,*,ザゾンビ,ザゾンビ

MeCabのドキュメントに詳しく書いてありますが、mecab-dict-index -d "C:\Program Files (x86)\MeCab\dic\ipadic" -u C:\path\to\animation.dic animation.csv して、res <- RMeCabC(text, dic="C:/path/to/animation.dic")と言うように使います。

トーラス世界の感染進行モデル

GitHubを使えって絶対に言うな!!!

えーっと、BCGワクチン接種が新型コロナウイルス感染症(COVID-19)に有効説検証論文に付随した話で、同一世界にワクチン接種ありとなしの集団がいたらどうなんの? — と言う御題があって、(ドラクエのマップのような縦横がつながっている)トーラス世界にマス目を区切って、そこに接種ありと接種なしの人々をランダムに配置して、隣接したマス目にいる人からのみ影響を受けるようなトイモデルで、どうなるか検討してみました。

大事なポイントなので強調しますが、以下のような感じで、乱雑にワクチン有りと無しを詰め込みます。

まぁ、ワクチンに効果があればワクチンを打つ方が感染しないですし、集団免疫獲得につれて効果量が特に落ちると言うわけでもないです。ただし、ワクチン有効率は最初が大きく中盤で落ち着きます。


しかし、詳しく特異なモデルをつくってもSIRモデルに近づくんですよね。やはり古典的な微分方程式モデルは強い。

set.seed(202004)

# 感染率β
beta_0 <- 0.50 # ワクチン非接種者
beta_1 <- 0.25 # ワクチン接種者

# 回復率γ
gamma <- 0.19

# ワクチン接種者の比率
theta <- 0.5

# トーラス状の空間を表すn行n列の行列を複数つくる
n <- 100

# 計算する期間
t <- 250

# ワクチン接種者は1で、他は0とする
shot <- numeric(n^2)
shot[order(runif(n^2))[1:as.integer(n^2*theta)]] <- 1

# ワクチン接種の有り無しの人数を数えておく
obs <- tapply(shot, shot, length)

# ワクチン接種の有無を表す行列
m_shot <- matrix(shot, n)

# 最初に感染しているのは1人
inf <- numeric(n^2)
inf[order(runif(n^2))[1]] <- 1

# 感染者を表す行列
m_inf <- matrix(inf, n)

# 暴露レベルを表す行列を返す
m_exposed <- function(m_inf){
  # 感染者の左側を1とする行列(トーラスなのに注意)
  inf <- c(m_inf)
  m_l <- matrix(c(inf[-(1:n)], inf[1:n]), n)
  
  # 感染者の右側を1とする行列
  inf <- c(m_inf)
  m_r <- matrix(c(inf[((n^2-n+1): n^2)], inf[-((n^2-n+1): n^2)]), n)
  
  # 感染者の上を1とする行列
  inf <- c(t(m_inf))
  m_u <- matrix(c(inf[-(1:n)], inf[1:n]), n, byrow=TRUE)
  
  # 感染者の下を1とする行列
  inf <- c(t(m_inf))
  m_b <- matrix(c(inf[((n^2-n+1): n^2)], inf[-((n^2-n+1): n^2)]), n, byrow=TRUE)

  # 0~4の値が返る
  m_l + m_r + m_u + m_b
}

# 完全な免疫を獲得した回復者を表す行列
m_rcv <- matrix(numeric(n^2), n)

# データフレームを作成する
df1 <- data.frame(
  t = 1:t,
  pct_infected_shot = numeric(t),
  pct_infected_no_shot = numeric(t),
  number_of_infected = numeric(t),
  VE = numeric(t)
)

for(i in 1:t){
  # 感染確率を表す行列
  # = ウイルスに暴露されている×社会的距離×{世代に応じたβ}×非回復者
  p <- (beta_0 - c(m_shot) * (beta_0 - beta_1)) * abs(c(m_rcv) - 1)
  # 周囲にいる感染者の人数に応じて、感染確率が上がる
  m_prob <- matrix(1 - (1 - p)^c(m_exposed(m_inf)), n)

  # 感染を判定し、感染者行列を更新
  m_inf[runif(n^2) < c(m_prob)] <- 1
  
  # 回復判定をし、感染者行列と免疫保持者行列を更新
  flag <- (runif(n^2) < gamma)&(m_inf==1)
  m_inf[flag] <- 0
  m_rcv[flag] <- 1

  # 感染経験者は、現在感染者もしくは回復者
  infected <- c(m_inf)|c(m_rcv)
  # ワクチン接種者と非接種者ごとに感染経験者を集計
  noi <- tapply(infected, c(m_shot), sum)
  # ワクチン接種者と非接種者ごとに感染経験率を出す
  pct <- noi/obs
  df1$number_of_infected[i] <- sum(noi)
  df1$pct_infected_no_shot[i] <- pct[1]
  df1$pct_infected_shot[i] <- pct[2]
  df1$VE[i] <- 1 - pct[2]/pct[1]
}

title <- sprintf("トーラス世界の感染進行モデル(%d行%d列,回復率%.2f,接種率%.2f)", n, n, gamma, theta)

windowsFonts(Meiryo = windowsFont("Meiryo"))
par(family="Meiryo", mar=c(5, 4, 4, 2), mgp=c(2, 1, 0), bg="white")

plot(df1$pct_infected_no_shot, col="red", type="l", lwd=2, axes=FALSE, xlab="", ylab="", lty=2, ylim=c(0, 1), main=title)
lines(df1$pct_infected_shot, col="blue", lwd=2)

y_at <- seq(0, 1, 0.25)
y_labels <- sprintf("%d%%", 100*y_at)
y_space <- abs(par()$usr[1] - par()$usr[2])/25
axis(2, at=y_at, labels=y_labels, pos=- y_space, las=2)

x_at <- seq(1, t, length.out=10)
x_labels <- sprintf("%.0f期", x_at)
x_space <- abs(par()$usr[3] - par()$usr[4])/40
axis(1, at=x_at, labels=x_labels, pos=-x_space, las=2)

legend("bottomright", lty=c(2, 1), lwd=c(2, 2), pch=c(-1, -1), col=c("red", "blue"), legend=c(sprintf("ワクチン接種無し(感染率β=%.2f)", beta_0), sprintf("ワクチン接種有り(感染率β=%.2f)", beta_1)), bty="n", y.intersp=1.25)

dev.copy(png, sprintf("%sβ0=%.2f,β1=%.2f.png", title, beta_0, beta_1), width=640, height=480, type="cairo", bg="white")
dev.off()


t2 <- as.integer(t/3)
df2 <- df1[1:t2,]
title <- sprintf("トーラス世界の感染進行モデル(%d行%d列,回復率%.2f,接種率%.2f)", n, n, gamma, theta)

windowsFonts(Meiryo = windowsFont("Meiryo"))
par(family="Meiryo", mar=c(5, 4, 4, 2), mgp=c(2, 1, 0), bg="white")

ylim <- c(0, ceiling(max(df2$pct_infected_no_shot)*10)/10)

plot(df2$pct_infected_no_shot, col="red", type="l", lwd=2, axes=FALSE, xlab="", ylab="", lty=2, ylim=ylim, main=title)
lines(df2$pct_infected_shot, col="blue", lwd=2)

y_at <- seq(0, ylim[2], length.out=5)
y_labels <- sprintf("%.0f%%", 100*y_at)
y_space <- abs(par()$usr[1] - par()$usr[2])/25
axis(2, at=y_at, labels=y_labels, pos=- y_space, las=2)

x_at <- seq(1, t2, length.out=10)
x_labels <- sprintf("%.0f期", x_at)
x_space <- abs(par()$usr[3] - par()$usr[4])/40
axis(1, at=x_at, labels=x_labels, pos=-x_space, las=2)

legend("topleft", lty=c(2, 1), lwd=c(2, 2), pch=c(-1, -1), col=c("red", "blue"), legend=c(sprintf("ワクチン接種無し(感染率β=%.2f)", beta_0), sprintf("ワクチン接種有り(感染率β=%.2f)", beta_1)), bty="n", y.intersp=1.25)

dev.copy(png, sprintf("%sβ0=%.2f,β1=%.2f,前半.png", title, beta_0, beta_1), width=640, height=480, type="cairo", bg="white")
dev.off()

Rでガブリエル比較区間を計算する方法

データを男女などで二つ以上の群にわけて、比較したいときは色々とあります。平均や分散の違いをt検定やF検定などで比較するのが教科書的かつ堅実な方法ですが、統計学に必ずしも詳しくない人々に説明する場合は、ビジュアライゼーションも考えていかないといけません。

この目的でよく使われるのは信頼区間によるエラーバーがついたグラフですが、実は往々にして誤解されています*1。2つの群のエラー・バーは重なっているのに、t検定で統計的に有意な差がでるときがあり、群を複数に分割したときの多重比較のための補正も行われません。比較対象が群ではなくて0などある一点であれば、信頼区間の厳密な説明*2とはちょっと違う使い方ではあるものの、t検定と合致した結果が得られるわけですが。

こういうわけで、群と群をグラフで比較するのには、信頼区間は使わない方が望ましいです。ガブリエル比較区間Gabriel (1978))と言う手法があるので、こちらを使うことをお勧めします。多重比較のための補正も行ってくれます。rgabrielと言うパッケージがあるので、使ってみましょう。

# 信頼区間のエラーバーが重なり、ガブリエル比較区間ではそうでないデータセット
df1 <- data.frame(category = c("C","C","C","C","C","C","T-1","T-1","T-1","T-1","T-1","T-1"), value = c(0.117,-0.718,0.7199,-0.2422,-0.3391,-0.2223,1.0094,1.4306,0.099,1.2839,0.403,0.5616))

# 有意水準(片側)
alpha <- 0.025

# ライブラリを呼び出す.インストールしていない場合は install.packages("rgabriel")を
library(rgabriel)

# 信頼区間とガブリエル比較区間を計算
attach(df1)
df2 <- data.frame(
  category=levels(category),
  count = as.numeric(tapply(value, category, length)),
  mean = as.numeric(tapply(value, category, mean)),
  sd = as.numeric(tapply(value, category, sd))
)
df2$se <- df2$sd/sqrt(df2$count-1)
df2$t_upper <- df2$mean + df2$se*qt(1 - alpha, df=df2$count-1)
df2$t_lower <- df2$mean - df2$se*qt(1 - alpha, df=df2$count-1)
df2$g_upper <- df2$mean + rgabriel(df1$value, df1$category, alpha)
df2$g_lower <- df2$mean - rgabriel(df1$value, df1$category, alpha)
detach(df1)

データフレームdf2の中のt_upperとt_lowerが信頼区間、g_upperとg_lowerがガブリエル比較区間になります。信頼区間では群Cのt_upperは群T-1のt_lowerよりも大きいですが、g_upperはg_lowerよりも小さくなっています。グラフを描くと、以下のようにガブリエル比較区間の方が縮まっているのが分かります。

f:id:uncorrelated:20190409102836p:plain
ガブリエル比較区間と信頼区間

t検定をかけて、どちらが正しいかを見てみましょう。

t.test(df1$value[df1$category=="C"], df1$value[df1$category=="T-1"])

両側検定で5%有意で帰無仮説が棄却され、群Cと群T-1に差があることが分かります。なお、多重比較のための補正は手法によって結果が変わるので、検定結果と整合的になるとは限らないことには注意してください。

グラフのソースコード

gabriel.plotがあるので、plotすればグラフになるのですがビジュアル的にイマイチだったので、以下のように描画しました。

windowsFonts(Meiryo = windowsFont("Meiryo UI"))
par(family="Meiryo", mar=c(2, 1, 1, 1), mgp=c(2, 1, 0))

unit <- 0.1
margin <- 0.05 # 信頼区間のエラーバーとガブリエル比較区間のエラーバーの隙間の幅
ylim <- c(unit*floor(min(df2$t_lower, df2$g_lower)/unit), unit*ceiling(max(df2$t_upper, df2$g_upper)/unit))
xlim <- c(c(0.5, length(df2$mean) + 0.5))
plot((1:length(df2$mean)) - margin/2, df2$mean, type="p", ylim=ylim, xlim=xlim, axes=FALSE, xlab="", ylab="", pch=19, col="red")
arrows(1:length(df2$mean) - margin/2, df2$g_lower, 1:length(df2$mean) - margin/2, df2$g_upper, code=0, col="red", lwd=2)
arrows(1:length(df2$mean) + margin, df2$t_lower, 1:length(df2$mean) + margin, df2$t_upper, code=0, lwd=2, lty=2)
points(1:length(df2$mean) + margin, df2$mean, pch=19)

abline(h=df2$t_upper[1], lty=2)
abline(h=df2$g_upper[1], col="red", lty=3)

abline(h=df2$t_lower[1], lty=2)
abline(h=df2$g_lower[1], col="red", lty=3)

axis(1, at=c(xlim[1], 1:length(df2$mean), xlim[2]), labels=c("",as.character(df2$category),""), pos=ylim[1], las=2)

legend("topleft", c(sprintf("%.0f%%信頼区間", (1-2*alpha)*100), "ガブリエル比較区間"), col=c("black", "red"), lty=c(2, 1), lwd=c(2, 2), bg='white', box.col = 'black', bty="n")

メタアナリシス

rgabriel関数のソースコードを見ていると、群を表すファクターごとの標本標準偏差のベクトル(s)と、ファクターごとのサブサンプル・サイズのベクトル(n)と、有意水準(a)があればガブリエル比較区間は計算できるので、複数の分析のデータをつなぐことも簡単にできると思います。

*1:ダメな統計学」第6章

*2: 真の平均値が観測値だとしたときに所定の観測数のデータをとると、観測される平均値が信頼区間に収まる確率が…と言うようなものです。

クロス集計表のセルごとの有意確率を出してみよう

他人への質問を眺めていても学ぶことはあります。昨日は「クロス集計の独立性検定を示した表で、セルごとにアスタリスクがあるのは何故か?」と言う質問がされていましたが、私も理由が思いつきませんでした。以下のような頻度表、時折、見てきた気がするのですが、どうやって個々のセルの有意性を出しているのでしょうか?

JPN USA CHN
MEAT 10 13 6
FISH 9 18** 3**
BEAN 20 15** 31**

Note) ** and * are statistically significant at the 1% and 5% level, respectively.

独立性検定(全体)

まずはχ二乗検定をしてみましょう*1
頻度が入った行列mを作ってchisq.test(m)で終わるのですが、今回は教科書的にチマチマと計算してみます。ただし、1%、5%棄却粋のχ二乗値を出して比較することはせずに、P値は計算してもらいます。

# 頻度の行列をつくる
m <- matrix(c(10,9,20,13,18,15,6,3,31), 3)

# 見やすいように行名と列名をつけておく
colnames(m) <- c("JPN", "USA", "CHN")
rownames(m) <- c("MEAT", "FISH", "BEAN")

# 列ごとの合計
apply(m, 2, sum)

# 行ごとの合計
apply(m, 1, sum)

# 全体の合計
sum(m)

# 自由度を計算する
degf <- (ncol(m) - 1)*(nrow(m) - 1)

# 行(MEAT, FISH, BEAN)ごとにの期待確率
exp_p <- apply(m, 1, sum)/sum(m)

# 列(JPN, USA, CHN)ごとに期待確率exp_pをかけて、期待値の行列をつくる
exp_m <- sapply(apply(m, 2, sum), function(x){
  x * exp_p
})

# 観測値の行列と期待値の行列のセルごとの差分の二乗を合計して期待値で割り、χ二乗値を求める
chisq_v <- sum((m - exp_m)^2/exp_m)

# P値を求める
pchisq(chisq_v, degf, lower.tail=FALSE)

当然ですが、0.0008233372と、一つしかP値は出ません。

セルごとに検定

残差の二乗値を期待値で割ったものの合計がχ二乗分布に従うことから検定をかけたわけですが、今度は個々のセルの調整済み*2残差が正規分布に従うことから検定をかけます。
これも2*pnorm(abs(chisq.test(m)$stdres), lower.tail=FALSE)で終わるのですが、チマチマとやっていきましょう。

# (1-列合計/合計)(1-行合計/合計)の修正用係数の行列
adj_m <- matrix(rep(1-apply(m, 2, sum)/sum(m), ncol(m)) * rep(1-apply(m, 1, sum)/sum(m), each=nrow(m)), nrow(m), byrow=TRUE)

# 標準化残差
stdres <- (m - exp_m)/sqrt(c(adj_m) * c(exp_m))

# セルごとのP値をまとめた行列をつくる
m_p_value <- 2*pnorm(abs(stdres), lower.tail=FALSE)

今回の結果は以下のようになりました。

JPN USA CHN
MEAT 0.6633 0.3064 0.1362
FISH 0.8707 0.0025 0.0030
BEAN 0.8189 0.0006 0.0001

多重比較のためのP値の調整

期待値1%の事象でも、100回繰り返すと64%の確率で1度は生じることになります。この理屈でクロス集計表が細かく切られている場合、偽陽性が出やすくなります。多重比較のために、P値を調整しましょう。ここはライブラリにお任せです。

m_adj_p_value <- matrix(p.adjust(m_p_value), nrow(m_p_value))
dimnames(m_adj_p_value) <- dimnames(m) # 見やすいように行名と列名をコピーしておく

p.adjustでとれる調整法は多彩ですが、デフォルトはHolm法になっています。

JPN USA CHN
MEAT 1 1.0000 0.6812
FISH 1 0.0176 0.0183
BEAN 1 0.0045 0.0013

今回は解釈を左右するような変化はないですが、P値が大きくなっていることに注意してください。FISH行の1%有意が5%有意になりました。なお、この手のP値の調整は偽陰性が大きくなる傾向がある*3そうです。

参考文献

本稿の作成において「クロス集計表の有意差検定 | 象と散歩」と「Rでクロス集計表の残差分析 - bob3’s blog」を参考にしました。

*1: サンプルサイズが小さく、セルの度数が一つでも5未満の場合は、多項分布を前提にしたFisherの正確検定を使うべきと言う話もありますが、ここでは無視します。

*2: 自由度を調整する作業になると思いますが、Agresti (2018) などで確認してください。Rのヘルプの雰囲気ではたぶん、理由が分かるはずです。たぶん(´・ω・`)ショボーン

*3: Bonferroni法、Holm法、False Discovery Rate | 大阪大学腎臓内科

行列を使った計算で、SURの感じを掴んでみよう

計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ないSUR *1 をふと思い出したので、手順を確認してみました。

1. 2つの方程式からなるSURモデル

理屈はGreeneのEconometric Analysisに詳しく載っているので参照して欲しいのですが、SURの概要を説明します。
以下の連立方程式を同時推定することを考えます*2

{\displaystyle 
\begin{eqnarray}
  \left\{
    \begin{array}{l}
		y_{1} = X_{1}\beta_1 + \epsilon_{1} \\
		y_{2} = X_{2}\beta_2 + \epsilon_{2}
    \end{array}
  \right.
\end{eqnarray}
}

 yは被説明変数のベクトル、 Xは説明変数の行列、 \betaは係数、 \epsilonは誤差項のベクトルです。添字の数字は方程式を表します。
一つ一つ推定しても良い気がしますが、 \mbox{COV}(\epsilon_1, \epsilon_2) \neq 0の場合は推定量を改善する余地があります。

例えば、天候データが無いときに2つの品種の収穫高と肥料の関係を推定するとして、日照量や気温など2つの方程式に同時に同じような影響を与える要因がある場合、 \mbox{COV}(\epsilon_1, \epsilon_2)という情報も使う方が推定結果は真の値に近づきます。

さて、最初は \mbox{COV}(\epsilon_1, \epsilon_2)の情報が無いので、以下をOLSで推定します。

{\displaystyle 
	\begin{pmatrix}
		y_{1} \\
		y_{2}
	\end{pmatrix}
	=
	\begin{pmatrix}
		X_{1} & 0 \\
		0 & X_{2}
	\end{pmatrix}
	\begin{pmatrix}
		\beta_{1} \\
		\beta_{2}
	\end{pmatrix}
	+
	\begin{pmatrix}
		\epsilon_1 \\
		\epsilon_2
	\end{pmatrix}
}

これで、以下の分散共分散行列 \Sigmaがつくれるようになりました。

{\displaystyle 
	\Sigma= 
	\begin{pmatrix}
		\mbox{COV}(\epsilon_1, \epsilon_1) & \mbox{COV}(\epsilon_1, \epsilon_2) \\
		\mbox{COV}(\epsilon_1, \epsilon_2) & \mbox{COV}(\epsilon_2, \epsilon_2)
	\end{pmatrix}
}

これをウェイトに使って、不均一分散の補正をした回帰(FGLS)をかけます。
つまり、以下のように表記を簡素化して、

{\displaystyle 
	y = X\beta + \epsilon
}

以下のようにGLS推定量を求めます。


\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
{\displaystyle 
	\hat{\beta} = \argmin_{\beta} \space  (y - X \beta)^t ( \Sigma^{-1} \otimes I_n ) (y - X \beta)
}

ここで、 \otimesクロネッカー積、 I_nはn行n列の対角行列、 \Sigma \otimes I_nは共分散を対角成分にとる細胞で構成される行列です。 nは観測数です。後述するコマンドを見れば、何が起きるかは一目瞭然だと思います。

この式を満たすために、 \beta微分して一階条件を出して整理すると、

{\displaystyle 
	\hat{\beta} = (X^t ( \Sigma^{-1} \otimes I_n ) X)^{-1} X^t ( \Sigma^{-1} \otimes I_n ) y
}

となります。 (\Sigma \otimes I_n)^{-1} = \Sigma^{-1} \otimes I_nに注意してください。

2. RによるSUR推定

systemfitパッケージを使えば良いのですが、前節で説明したモデルとの対応関係が分からなくなるので、行列を使って計算します。

#
# データセットの作成
#
set.seed(20181103)
n <- 50 # 標本サイズが小さい方が、OLSとSURの差が出る

l <- rnorm(n, sd=2) # 共分散を持つようにする
e1 <- rnorm(n, sd=1) + l
e2 <- rnorm(n, sd=1) + l

x1 <- runif(n, min=0, max=2)
x2 <- runif(n, min=-2, max=2)
y1 <- 1 + 2*x1 + e1
y2 <- 3 + 2*x2 + e2

X1 <- matrix(c(rep(1, n), x1, rep(0, 2*n)), n, 4)
X2 <- matrix(c(rep(0, 2*n), rep(1, n), x2), n, 4)
X <- rbind(X1, X2)
y <- c(y1, y2)

#
# 以下、実際のSUR推定量の計算
#
# 1段階目のOLS推定量を計算
beta_ols <- solve(t(X)%*%X)%*%t(X)%*%y

# 誤差項を計算
r_e <- y - X %*% beta_ols
r_e1 <- r_e[1:n]
r_e2 <- r_e[(n+1):(2*n)]

# ∑を計算
sigma <- matrix(c(cov(r_e1,r_e1), cov(r_e1,r_e2), cov(r_e1,r_e2), cov(r_e2,r_e2)), 2, 2)

# SUR推定量を計算
beta_sur <- solve(t(X) %*% (solve(sigma) %x% diag(n)) %*% X) %*% t(X) %*% (solve(sigma) %x% diag(n)) %*% y

なお、beta_ols - beta_surで推定量にどれぐらい違いがあるか分かるので、nの値を変えて試してみてください。t値やP値はOLSと同様に計算できます。また、sigma %x% diag(n)とすると、 n  \times nの巨大行列ですが、 \Sigma \otimes I_nがどうなっているか分かります。

*1: すぐに同時性ガガガ…と言われるので、GMMか3SLS、もしくは全情報最尤法を使うのが一般的で、SUR自体は3SLSを学ぶ前の予備知識的な扱いが多いと思います。

*2:もっとも単純なモデルを考えましたが、二つの式で、係数が同じ値を持つような制約を置くこともできます。

相関行列から一定以上の大きさの変数の組み合わせを抽出

100や1000もある説明変数同士の相関係数を見るはだるいと言う話があって、確かに100×100や1000×1000どころか、10×10の行列でも目視で確認すると見落としが出そうです。相関係数0.5以上をリストするコードを書いてみましょう。

# テスト用のデータフレームを作る
set.seed(20180925)

n <- 100 # 変数の長さ
m <- 100 # 変数の種類
df1 <- data.frame(x001 = rnorm(n))
for(i in 2:m){
  vn <- sprintf("x%03d", i)
  df1[vn] <- rnorm(n)
  # 以下、ランダムに、既に生成済みの変数からの影響を加える
  xi <- as.integer(runif(1, max=i))
  if(0<xi){
    xn <- sprintf("x%03d", xi)
    beta <- runif(1, min=-1, max=1)
    df1[vn] <- df1[vn] + beta*df1[xn]
  }
}

# だるさを確認
cor(df1)

# 機械的相関係数が0.5より大なものを選ぶ
t <- 0.5
n <- names(df1)
l <- length(n)
s <- 1:l^2-1
flag <- c(abs(cor(df1))>t) # 相関係数がtより大なものはTRUE
flag[seq(1, l^2, l + 1)] <- FALSE # 自己相関は排除

# 結果表示
pair <- sprintf("cor(%s,%s)>%.1f", n[floor(s/l) + 1], n[s %% l + 1], t)
pair[flag]

9900の組み合わせを目視するよりはましですが、116個も出てくるとどれを落とそうか迷いそうですね。こういう状態で回帰分析をするときは、L1/L2正則化項を加えるか、主成分分析などで変数をまとめた方が良さそうです。やはり見る意味はないかも知れません。

よくある操作変数法(IV)のための練習データとOLSによる推定バイアス

内生性について上手くツイートでは問題点を説明できないので、念頭に置いているコードを出します。

#
# 教科書例の需要と供給のデータを作る
# 真のモデル:
# S = 2 + 3*p + 4*z + ν
# D = 1 - 1*p + μ
# S = D
# p:価格, S:供給, D:需要, z:気候か何か, νとμ:誤差項
# → 価格について式を整理してみると、D=Sと言う式から、P = (2 - 1)/(-1 - 3) + 4・z/(-1 - 3) + (ν - μ)/(-1 - 3)になって、Pに推定する需要関数の誤差項μが含まれていることになり、OLSの条件を満たしません。
#
set.seed(20180411)
obs <- 300 # 観測数(増やしてバイアスの変化をチェック!)
a0 <- 1
a1 <- -1
b0 <- 2
b1 <- 3
b2 <- 4
z <- runif(obs, min=0, max=3)
mu <- rnorm(obs, mean=0, sd=2)
nu <- rnorm(obs, mean=0, sd=1)
p <- (b0 - a0 + b2*z + nu - mu)/(a1 - b1)
df1 <- data.frame(
  p = p,
  d = a0 + a1*p + mu,
  s = b0 + b1*p + b2*z + nu
)
rm(p)

#
# OLSで推定してみる
# → バイアスが入る
#
r_ols <- lm(d ~ p, data=df1)


#
# クロスセクションチェック
# → バイアスは同じまま
#
r_ols_css1 <- lm(d ~ p, data=df1[seq(1, obs, 2), ])
r_ols_css2 <- lm(d ~ p, data=df1[seq(2, obs, 2), ])

#
# tslsパッケージを使えば済む操作変数法
# → バイアスは入らない
#
attach(df1)
zm <- matrix(c(rep(1, obs), c(z)), obs, 2)
xm <- matrix(c(rep(1, obs), c(p)), obs ,2)
iv_estimated_a <- solve(t(zm) %*% xm) %*% (t(zm) %*% d)
detach(df1)

#
# 推定結果の中身を比較してみる
#
print_dc <- function(label, a){
  sprintf("%s: D = %.04f%s%.04fp", label, a[1], ifelse(a[2]>=0,"+",""), a[2])
}

paste(c(
  print_dc("OLS", coef(r_ols)),
  print_dc("CSS01", coef(r_ols_css1)),
  print_dc("CSS02", coef(r_ols_css2)),
  print_dc("IV", iv_estimated_a)))