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

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

	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そうです。

参考文献

本稿の作成において「クロス集計表の有意差検定 | 象と散歩」と「Ｒでクロス集計表の残差分析 - 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。

は被説明変数のベクトル、は説明変数の行列、は係数、は誤差項のベクトルです。添字の数字は方程式を表します。
一つ一つ推定しても良い気がしますが、の場合は推定量を改善する余地があります。

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

さて、最初はの情報が無いので、以下をOLSで推定します。

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

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

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

ここで、はクロネッカー積、はn行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)]

# ∑を計算
df <- n - 2 # 自由度は複数の方程式の最も小さい値になる
sigma <- matrix(c(sum(r_e1^2)/df, sum(r_e1*r_e2)/df, sum(r_e1*r_e2)/df, sum(r_e2^2)/df), 2, 2)

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

# 分散共分散行列を計算
r <- (y - X %*% beta_sur)
vcov <- solve(t(X) %*% (solve(sigma) %x% diag(n)) %*% X) # Usual Variance Matrix (Wooldridge (2002) p.161)

なお、beta_ols - beta_surで推定量にどれぐらい違いがあるか分かるので、nの値を変えて試してみてください。t値やP値はOLSと同様に計算できます。また、sigma %x% diag(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)))

推定結果は

OLS:
CSS01:
CSS02:
IV:

となり、OLSは真のモデルとは大きく異なる推定結果をだし、クロスチェックもそれを検知しない一方、操作変数法が上手く機能していることが分かります。

FE-IV練習用データ生成から、one-wayクラスター頑強標準誤差の計算まで

dとsの式でpが内生変数、zが操作変数、iが個体を表す番号ですが、こんなんで。

noi <- 20 # 個体数
t <- 5 # 観測期間
obs <- t*noi
i <- rep(1:(noi), each=t)
fe <- runif(noi, min=0, max=100)
a0 <- rep(fe, each=t)
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)
d <- a0 + a1*p + mu
s <- b0 + b1*p + b2*z + nu

within変換をかけてIVで推定

within_transfer <- function(x, i){
  m <- tapply(x, i, mean)
  x - rep(m, each=t)
}

w_d <- within_transfer(d, i)
w_z <- within_transfer(z, i)
w_p <- within_transfer(p, i)

zm <- matrix(c(c(w_z)), obs, 1)
xm <- matrix(c(c(w_p)), obs ,1)
estimated_a1 <- solve(t(zm) %*% xm) %*% (t(zm) %*% w_d)

観測数obsが2000ぐらい無いと誤差多し

Stata風に切片項を計算する

estimated_mu <- w_d - estimated_a1 %*% w_p

ma_d <- mean(d)
mi_d <- rep(tapply(d, i, mean), each=t)

ma_z <- mean(z)
mi_z <- rep(tapply(z, i, mean), each=t)

ma_p <- mean(p)
mi_p <- rep(tapply(p, i, mean), each=t)

estimated_a0 <- (d - mi_d + ma_d) - (p - mi_p + ma_p)*c(estimated_a1) - c(estimated_mu)

切片項のP値が欲しい場合

StataのFAQを見ると、within変換をかけた変数に全体の平均値を加算してから推定をすることで、切片項の有意性を出していました。

w_d <- within_transfer(d, i) + mean(d)
w_z <- within_transfer(z, i) + mean(z)
w_p <- within_transfer(p, i) + mean(p)

zm <- matrix(c(rep(1, obs), c(w_z)), obs, 2)
xm <- matrix(c(rep(1, obs), c(w_p)), obs ,2)
estimated_a <- solve(t(zm) %*% xm) %*% (t(zm) %*% w_d)

分散共分散行列と標準誤差の計算は以下です。自由度の計算でパネルの数（個体種類）も加味しないと、within推定とdummy variable estimatorの標準誤差が一致しないことに注意しましょう。

df <- (obs - 2) - (noi - 1)
ssr <- sum((w_d - xm %*% estimated_a)^2)
s2 <- ssr/df
vcov <- s2*solve( t(xm) %*% zm %*% solve(t(zm) %*% zm) %*% t(zm) %*% xm )
se <- sqrt(diag(vcov))

one-wayクラスター頑強標準誤差の計算

気づくと一般化したone-wayのロバスト標準誤差（ここでは個体ごとの不均一分散を調整）の場合は以下です。「切片項のP値が欲しい場合」の続きとして実行できます。

# 時点を示す番号（two-way clustering用）
# ti <- rep(1:t, noi)

# 残差からウェイトΩを計算
residuals <- (w_d - xm %*% estimated_a)
omega <- matrix(0, obs, obs)
# Ωは対角成分としてt×tの細胞をnoi個とる
for(j in 1:obs){
  omega[,j] <- residuals[j]*residuals*(i[j]==i)
  # two-way clusteringのときは以下にする
  # omega[,j] <- residuals[j]*residuals*(i[j]==i|ti[j]==ti)
}

# plmパッケージのvcovHC(model, type = 'HC0') と同じ値
vcov_hc0 <- solve( t(xm) %*% zm ) %*% (t(zm) %*% omega %*% zm) %*% solve( t(zm) %*% xm )

# 自由度を調整
df <- nrow(zm) - ncol(zm)
dfcw <- df / (df - (noi - 1))
dfc <- (noi / (noi - 1))*((obs - 1)/df)
vcov <- dfc*vcov_hc0*dfcw

# 標準偏差を計算
se <- sqrt(diag(vcov))

IVとFGLSのあわせ技になっていますが、xmとzmを同じにすればIVでないのと同じになります。
なお、太田 (2013)と、パッケージを使わないようにコードは変えましたが、ストックホルム大学のMahmood Arai教授のレクチャーノートを参考にしました。

plmパッケージでのone-wayクラスター頑強標準誤差の計算の仕方

もっとも需要が大きそうな話を書き忘れていたので4年経過していますが追記します。Stata風に切片項を推定する方法は無さそうですが、自由度調整をすれば同様にクラスター頑強標準誤差になります。

library(plm)
# データフレームをつくる
df01 <- data.frame(i, t = rep(1:5, noi), d, s, p, z)
# パネルデータ分析用データフレームにする
pdf01 <- pdata.frame(df01, index=c("i", "t"))
# within推定をする
r_plm <- plm(d ~ p | z, data = pdf01, model = "within")

# 自由度を調整
# df <- with(r_plm, nrow(model)- (length(coefficients)) + 1)
# dfcw <- df / (df - (noi - 1)) # noi:観測数
# dfc <- (noi / (noi - 1))*((obs - 1)/df) # obs: サンプルサイズ
# vcov <- dfc*vcovHC(r_plm, type = 'HC0')*dfcw
# 
# 補正VCOVで標準誤差を計算
# summary(r_plm, vcov = vcov)
# 
summary(r_plm, vcov = vcovHC(r_plm, type = 'sss')) # これでStata互換になるそうです：https://blog.theleapjournal.org/2016/06/sophisticated-clustered-standard-errors.html

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

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

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

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))

分散拡大係数（VIF）は1.02弱と小さいです。

1/(1-summary(lm(x1 ~ x2))$r.squared)

多重共線性があるケース

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

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)

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

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

VIFは2.81強と大きくは無いですが、観測数が小さいため影響しています。

1/(1-summary(lm(x1 ~ x2))$r.squared)

多重共線性が深刻なケース

さすがにレアなのですが、真のモデルがのとき、と相関する変数を加えてを推定すると、の有意性が消えて係数が逆転し、の係数が有意になるデータをつくります。誤ったモデルの方が、自由度調整済み重相関係数は改善することに注意してください。

set.seed(737)
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)
e <- rnorm(n, sd=10)

# 説明変数の高い相関，高いVIF
sprintf("cor:%f, VIF:%f", cor(x1, x2), 1/(1-summary(lm(x1 ~ x2))$r.squared))

# 真のモデル/データ生成モデル
y <- 1 + 2*x1 + e

# 真のモデルに近い係数が有意に得られる
summary(lm(y ~ x1))

# x1の係数が逆転し有意性が消え、x2の係数が有意になる
summary(lm(y ~ x1 + x2))

サンプルサイズが小さければ小さいほど、誤差項の分散が大きければ大きいほど、VIFが高ければ高いほど色々とおきます。例ではR²は低いですが、R²が0.6ぐらいあっても、VIFが極端に高いと多重共線性が問題になる可能性はあります*1。

LOOCV

多重共線性が深刻なケースでLOOCVをかけるとx2だけで回帰した方がよくなります。

library(boot)
df01 <- data.frame(y, x1, x2)
frml <- c("y ~ x1", "y ~ x2", "y ~ x1 + x2")
delta <- matrix(NA, length(frml), 2)
rownames(delta) <- frml
colnames(delta) <- c("raw", "adjusted")
for(i in 1:length(frml)){
    delta[i, ] <- cv.glm(df01, glm(formula(frml[i]), data=df01))$delta
}
print(delta)

*1:adj-R²が0.58でも、VIFが89で引き起こせるケースがありました。100回やって7回ぐらいなので偶然の範疇と言えば、範疇なのかもですが。