Rで推定する回帰モデル入門
Sho Kuroda
May 3, 2017
Rで行う回帰分析入門．もともとR入門などに書いていた内容をこちらに移しました．パッケージのサンプルコードをコピー&ペーストしただけの部分も数多くあります．
Standard/classical regression
はじめに「y = Xb」の形式で記述される標準的な回帰モデルの推定を行います．
余談ですが，回帰モデルの説明変数選択の際には，互いに独立ではない／相関の高い変数（多重共線性）に注意する必要があります．前者については，カテゴリカル変数で全カテゴリーについて {0, 1} インディケータをモデルに投入することはもちろん，割合のデータ（組成データ）で複数カテゴリの割合を同時に投入すること，さらには地理空間データでは複数の距離指標を同時に投入することなども該当します（cf. Wang et al. (2015, IJCESCAE)）．さらに因果構造を考慮しない変数の選択によって共線性が生じたり，あるいは媒介／中間変数の不適切な特定化によって因果効果の推定にバイアスが生じることもあります（cf. 星野 (2009, 岩波)）．
OLS estimation for linear model / 線形モデルのOLS推定
\[ y = X \beta + \varepsilon, \varepsilon \sim N (0, \sigma^2) \]
まずは最小二乗法（OLS）．推定を行う関数 lm (fitting linear models) は標準で用意されています（OLS推定を行う lsfit という関数もありますが，忘れてよい）．チルダ ~ の前後に被説明変数と説明変数をそれぞれ指定するだけでOK．
y1 <- c(1, 2, 4, 5, 7, 8) # 被説明変数（応答変数）にするデータを適当に用意
x1 <- c(3, 5, 8, 8, 11, 10) # 説明変数
plot(x = x1, y = y1) # 描画
lm1 <- lm(formula = y1 ~ x1) # OLS推定
lm1 # 回帰係数しか表示されない##
## Call:
## lm(formula = y1 ~ x1)
##
## Coefficients:
## (Intercept) x1
## -2.0110 0.8681summary(object = lm1) # 標準誤差や決定係数まで自動で計算してくれる##
## Call:
## lm(formula = y1 ~ x1)
##
## Residuals:
## 1 2 3 4 5 6
## 0.40659 -0.32967 -0.93407 0.06593 -0.53846 1.32967
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.0110 1.0609 -1.896 0.13090
## x1 0.8681 0.1328 6.538 0.00283 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8957 on 4 degrees of freedom
## Multiple R-squared: 0.9144, Adjusted R-squared: 0.893
## F-statistic: 42.75 on 1 and 4 DF, p-value: 0.002828余談ですが，実は「 formula = 」部分はわざわざ書かなくても大丈夫．算術平均を計算する関数 mean を用いるときは，わざわざ「 mean(x = hoge) 」のように「 x = 」を含めては書きませんでしたよね，それと同じです（もちろん指定してもよいですが，明示するまでもなくほとんど自明な場合は省略してよいでしょう）．いずれの引数（ひきすう）に対応するかが省略されている場合は，関数の定義されている順番に使用されます．summary 関数の「 object = 」も同様です（これ以降は省略します）．
閑話休題．説明変数として変数の二乗を使用したいときは？単に「 lm(formula = y1 ~ x1^2) 」とするだけではダメで，関数 I を噛ませる必要があります．lm2 <- lm(formula = y1 ~ I(x1^2))説明変数が複数ある場合（重回帰分析のケース）は，プラス記号 + で繋げばよい．
lm3 <- lm(formula = y1 ~ x1 + I(x1^2))Rに標準で用意されているデータセットでも試してみましょう．データセットに含まれる変数のうち被説明変数以外の全てを説明変数として使用する場合は，formula引数の説明変数指定箇所をピリオド . で省略することができます． swiss は1988年頃のスイスにおける出生率 fertility と社会経済指標を地域単位で整備した cross section データです．
data(swiss) # データを呼び出す呪文（実は今回のデータはわざわざ呼び出さなくても使えるのですが，一応）
plot(x = swiss$Agriculture, y = swiss$Fertility) # 男性の農業従事者割合と出生率の散布図
summary(lm(formula = Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data = swiss))##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education +
## Catholic + Infant.Mortality, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2743 -5.2617 0.5032 4.1198 15.3213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
## Agriculture -0.17211 0.07030 -2.448 0.01873 *
## Examination -0.25801 0.25388 -1.016 0.31546
## Education -0.87094 0.18303 -4.758 2.43e-05 ***
## Catholic 0.10412 0.03526 2.953 0.00519 **
## Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
## F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10summary(lm(formula = Fertility ~ . , data = swiss)) # 同上##
## Call:
## lm(formula = Fertility ~ ., data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2743 -5.2617 0.5032 4.1198 15.3213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
## Agriculture -0.17211 0.07030 -2.448 0.01873 *
## Examination -0.25801 0.25388 -1.016 0.31546
## Education -0.87094 0.18303 -4.758 2.43e-05 ***
## Catholic 0.10412 0.03526 2.953 0.00519 **
## Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
## F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10ここまで「データセット」などと不明瞭な表現をしてきましたが，上記 swiss は「data.frame」という型のデータです（オブジェクトの型は class 関数で確認可能）．data.frame型のデータの特定の変数（列）にアクセスするにはドル記号 $ を使用します（上記 plot 関数の中身を参照．行列と同様に swiss[, 1] のような列番号による指定も可能）．
さて，式が長くなってくると可読性を損ないますし，同じ式を何度も使用する際はコピー＆ペーストするのも非効率．実はこの部分もオブジェクトに代入できます（formula型）． update 関数によって回帰式の変数を加減できるのも便利です．
formula1 <- Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality
lm_swiss1 <- lm(formula = formula1, data = swiss)
lm_swiss2 <- update(object = lm_swiss1, formula. = . ~ . - Infant.Mortality) # 説明変数から "Infant.Mortality" を除去
summary(lm_swiss2)##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education +
## Catholic, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.7813 -6.3308 0.8113 5.7205 15.5569
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.05542 6.94881 13.104 < 2e-16 ***
## Agriculture -0.22065 0.07360 -2.998 0.00455 **
## Examination -0.26058 0.27411 -0.951 0.34722
## Education -0.96161 0.19455 -4.943 1.28e-05 ***
## Catholic 0.12442 0.03727 3.339 0.00177 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.736 on 42 degrees of freedom
## Multiple R-squared: 0.6498, Adjusted R-squared: 0.6164
## F-statistic: 19.48 on 4 and 42 DF, p-value: 3.95e-09formula2 <- update.formula(old = formula1, new = . ~ . - Catholic) # 回帰式のみ更新することも可能
summary(lm(formula = formula2, data = swiss))##
## Call:
## lm(formula = formula2, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.513 -5.011 -1.188 5.522 16.982
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.77314 11.62836 5.914 5.28e-07 ***
## Agriculture -0.12929 0.07485 -1.727 0.09144 .
## Examination -0.68799 0.22628 -3.040 0.00406 **
## Education -0.61965 0.17631 -3.515 0.00107 **
## Infant.Mortality 1.30710 0.40658 3.215 0.00251 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.796 on 42 degrees of freedom
## Multiple R-squared: 0.6444, Adjusted R-squared: 0.6105
## F-statistic: 19.02 on 4 and 42 DF, p-value: 5.413e-09ところでOLSEを行列・ベクトル表記で表せば beta = (X’ X)^(-1) X’ y です．列ベクトルの結合が cbind 関数，行列の転置が t 関数，逆行列が solve 関数，行列の積が %*% （* は要素同士の積）であることを思い出せば，以下のようにコーディングすることができます．
\[ \beta_{OLS} = (X' X)^{-1} X' y \]
y_swiss <- swiss$Fertility # 被説明変数
x_swiss <- cbind(1, swiss$Agriculture, swiss$Examination) # 説明変数（定数項の 1 を含む点に注目）
solve(t(x_swiss) %*% x_swiss) %*% t(x_swiss) %*% y_swiss # 回帰係数のみ## [,1]
## [1,] 94.60968827
## [2,] -0.09399751
## [3,] -1.19502875lm(Fertility ~ Agriculture + Examination, swiss) # 確認##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination, data = swiss)
##
## Coefficients:
## (Intercept) Agriculture Examination
## 94.610 -0.094 -1.195しかしこれでは回帰係数の標準誤差やモデルの決定係数は分かりません． summary(lm(…)) を実行した時のようにまとめて表示する関数を書いてみましょう．
\[ e = y - \hat{y}, \hat{\sigma^2} = \frac{RSS}{n-k}, V(\beta) = \hat{\sigma^2} (X' X)^{-1}, \frac{\hat{\beta}}{SE} \sim t, R^2 = 1 - \frac{RSS}{TSS} \]
ついでと言っては何ですが，誤差項の不均一分散に対して頑健（robust）な標準誤差も計算しておきましょう．パッケージに頼るなら car::hccm などが使用できます．
\[ V_{\rm{white}} (\beta) = (X' X)^{-1} X' \Omega X (X' X)^{-1}, \hat{\Omega} = \rm{diag} (e_1^2, \ldots, e_n^2) \frac{n}{n - k} \]
lm_ols <- function (y_ols, x_ols) {
n_ols <- nrow(x_ols) # 説明変数行列の行数（サンプルサイズ）
k_ols <- ncol(x_ols) # 説明変数行列の列数（定数項を含めた説明変数の数）
xxi <- solve(t(x_ols) %*% x_ols) # solve(crossprod(x = x_ols)) でもよい
beta_ols <- as.numeric(xxi %*% t(x_ols) %*% y_ols) # 回帰係数
residuals_ols <- y_ols - x_ols %*% beta_ols # 残差ベクトル
residuals_ols <- as.numeric(residuals_ols)
rss_ols <- sum(residuals_ols^2) # RSS
tss_ols <- sum((y_ols - mean(y_ols))^2) # TSS: 「var(y_ols) * (n_ols - 1)」でも同じ
S2_ols <- rss_ols / (n_ols - k_ols) # sigma^2 の不偏推定量
var_beta_ols <- S2_ols * solve(t(x_ols) %*% x_ols) # 回帰係数の推定量の分散共分散行列
se_ols <- sqrt(diag(var_beta_ols)) # 標準誤差
t_ols <- beta_ols / se_ols # t値
p_ols <- 2*(1 - pt(q = abs(t_ols), df = n_ols - k_ols)) # p値
# Whiteの不均一分散頑強（ロバスト）標準誤差
var_beta_white_hc0 <- xxi %*% t(x_ols) %*% diag(residuals_ols^2) %*% x_ols %*% xxi # original
var_beta_white_hc1 <- xxi %*% t(x_ols) %*% diag(residuals_ols^2) %*% x_ols %*% xxi * n_ols/(n_ols-k_ols)
se_white <- sqrt(diag(var_beta_white_hc1))
t_white <- beta_ols / se_white
p_white <- 2*(1 - pt(q = abs(t_white), df = n_ols - k_ols))
est_ols <- cbind(Est = beta_ols, SE = se_ols, t = t_ols, p = p_ols,
SE.white = se_white, t.white = t_white, p.white = p_white)
r2_ols <- 1 - (rss_ols / tss_ols) # 決定係数
r2_adj_ols <- 1 - (1 - r2_ols) * (n_ols - 1) / (n_ols - k_ols) # 自由度調整済み
list(sampleSize = n_ols, Estimates = est_ols, R2 = r2_ols, adjR2 = r2_adj_ols)
}計算してみる．
y_swiss <- swiss$Fertility # y
x_swiss <- cbind(1, swiss$Agriculture, swiss$Examination) # X
lm_ols(y_ols = y_swiss, x_ols = x_swiss)## $sampleSize
## [1] 47
##
## $Estimates
## Est SE t p SE.white t.white
## [1,] 94.60968827 7.82710709 12.087440 1.332268e-15 5.53351838 17.097565
## [2,] -0.09399751 0.08590323 -1.094226 2.798103e-01 0.06325394 -1.486034
## [3,] -1.19502875 0.24454695 -4.886705 1.400471e-05 0.19537608 -6.116556
## p.white
## [1,] 0.000000e+00
## [2,] 1.444008e-01
## [3,] 2.274423e-07
##
## $R2
## [1] 0.4326045
##
## $adjR2
## [1] 0.4068138summary(lm_swiss3 <- lm(Fertility ~ Agriculture + Examination, swiss)) # OLSE##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.4089 -6.3234 0.0577 6.3134 20.8937
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.6097 7.8271 12.087 1.41e-15 ***
## Agriculture -0.0940 0.0859 -1.094 0.28
## Examination -1.1950 0.2445 -4.887 1.40e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.621 on 44 degrees of freedom
## Multiple R-squared: 0.4326, Adjusted R-squared: 0.4068
## F-statistic: 16.77 on 2 and 44 DF, p-value: 3.85e-06sqrt(diag(car::hccm(model = lm_swiss3, type = "hc1"))) # ロバスト標準誤差## (Intercept) Agriculture Examination
## 5.53351838 0.06325394 0.19537608lm 関数のありがたみが分かります．
Multicollinearity / 多重共線性
まずはVIF (variance inflation factor; 分散拡大要因) の計算から．方法はいろいろあって…
- 群馬大・青木先生の tolerance 関数を使う → source(file = “http://aoki2.si.gunma-u.ac.jp/R/src/tolerance.R”, encoding = “euc-jp”)
- {car} パッケージの vif 関数を使う
- {DAAG} パッケージの vif 関数を使う
- {HH} パッケージの vif 関数を使う
ここでは {car} パッケージの vif 関数を試してみましょう．他のパッケージでも同名の関数が定義されているため，パッケージを同定する意味で car::vif のように書きます．
lm_swiss1 <- lm(Fertility ~ . , swiss)
car::vif(lm_swiss1)## Agriculture Examination Education Catholic
## 2.284129 3.675420 2.774943 1.937160
## Infant.Mortality
## 1.107542パッケージを使用しなくても，定義から以下のように計算できます．「 summary(lm(HOGE ~ . - Fertility, swiss))$r.squared 」は，説明変数のうち HOGE をそれ以外の説明変数に回帰した際の回帰モデルの決定係数です．
\[ \rm{VIF}_k = \frac{1}{1 - R_k^2} \]
1 / (1 - summary(lm(Agriculture ~ . - Fertility, swiss))$r.squared)## [1] 2.2841291 / (1 - summary(lm(Examination ~ . - Fertility, swiss))$r.squared)## [1] 3.67542# 以下同様上は「（定数項部分を除いた）説明変数行列 swiss[, 2:6] の相関係数行列 cor の逆行列 solve の対角成分 diag」（助詞の連続は見逃してください）と同値であることを思い出せば以下のように簡潔に計算できます（実は青木先生のコードもこの方法を使用している）．ただし完全な線形従属関係が生じている場合には計算できません．
diag(solve(cor(x = swiss[, 2:6])))## Agriculture Examination Education Catholic
## 2.284129 3.675420 2.774943 1.937160
## Infant.Mortality
## 1.107542OLSE with optim function / optim を用いた最小二乗法
最尤（ML）法によってモデルを推定する場合，汎用の最適化関数 optim を用いて対数尤度関数を最大化します．最小二乗推定量 OLSE は最適化関数など使わずとも解析的にパラメータが導出できますが，そうでない多くの場合においては optim を使わざるを得ません．そのための手慣らしとして，まずは最小二乗法を実装してみましょう．なお「汎用関数＝テンプレート関数」という意味ではこの用語を使用していません，念のため．
\[ \hat{\beta_{\rm{OLS}}} = \arg \min (e'e) = \arg \min \sum_i (y_i - x_i' \beta)^2 \]
optim 関数はパラメータの初期値 par と最小化する関数 fn を引数として取り，目的関数を最小にするパラメータ $par とそのときの目的関数の値 $value などを戻り値として返します．
ols <- function (y_ols, x_ols) { # 最小二乗法を実行する関数
# y_ols 被説明変数ベクトル
# x_ols 説明変数行列
min_rss <- function (param) { # param （ここでは回帰係数）の関数となっている点に注目
sum((y_ols - x_ols %*% param)^2)
}
k <- ncol(x_ols)
return(optim(par = rep(0, k), fn = min_rss)) # 「return」は省略できる
}
y_swiss <- swiss$Fertility # y（先程と同じ）
x_swiss <- cbind(1, swiss$Agriculture, swiss$Examination) # X
ols(y_ols = y_swiss, x_ols = x_swiss) # 戻り値のうち「$par」が回帰係数## $par
## [1] 94.63321661 -0.09419927 -1.19589201
##
## $value
## [1] 4072.741
##
## $counts
## function gradient
## 238 NA
##
## $convergence
## [1] 0
##
## $message
## NULL以下のように短く書くこともできますが，後々読み返すことを考慮すればそれほどお勧めできません．自分で書いたコードでも，数ヶ月経てばコメントなしでは読めなくなります．
optim(c(0, 0, 0), function (param) sum((y_swiss - x_swiss %*% param)^2))## $par
## [1] 94.63321661 -0.09419927 -1.19589201
##
## $value
## [1] 4072.741
##
## $counts
## function gradient
## 238 NA
##
## $convergence
## [1] 0
##
## $message
## NULLこれが分かると，最小絶対値法 LAD （beta = arg min |e|）も簡単に実装できます．
optim(c(0, 0, 0), function (param) sum(abs(y_swiss - x_swiss %*% param)))## $par
## [1] 99.1842110 -0.1833711 -1.1262065
##
## $value
## [1] 325.5182
##
## $counts
## function gradient
## 346 NA
##
## $convergence
## [1] 0
##
## $message
## NULLquantreg::rq(formula = Fertility ~ Agriculture + Examination, data = swiss) # 答え合わせ## Call:
## quantreg::rq(formula = Fertility ~ Agriculture + Examination,
## data = swiss)
##
## Coefficients:
## (Intercept) Agriculture Examination
## 99.1843305 -0.1833722 -1.1262108
##
## Degrees of freedom: 47 total; 44 residualMLE of linear model / 線形モデルの最尤推定
See also:
- 最尤法によるフィッティング （奥村先生）
回帰係数については OLSE も y ～ N(0, sigma^2) を仮定した最尤推定量 MLE も同じで解析的に解けますが，練習のために対数尤度関数の最大化によって導出してみましょう．
回帰係数から標準誤差，AICまでまとめて計算する関数を書きましょう．対数化せずに（尤度関数を prod(dnorm(x = Y - X %*% beta)) のように設定して）力ずくで解こうとしても，数値計算上の問題が生じて推定が難しくなります（まず無理）．ここで，最適化するパラメータ param は偏回帰係数 (beta0, beta1, …) だけではなく誤差分散 sigma^2 を含んだベクトルであることに注意しましょう．
\[ \ln L = l = - \frac{n}{2} \ln (2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - x_i' \beta)^2 \]
optim 関数はデフォルトでは指定された関数（今回は尤度関数）の最小化を行うため，control 部分で最大化するように修正する必要があります．
ml_gauss <- function (y_gauss, x_gauss) {
LL_gauss <- function (param) { # 対数尤度関数
- (n / 2) * log(2 * pi * param[k+1]) -
sum((y_gauss - x_gauss %*% param[1:k])^2) / (2 * param[k+1]) # param[k+1] は sigma^2 の推定量に相当
}
n <- nrow(x_gauss) # サンプルサイズ
k <- ncol(x_gauss) # 説明変数の数（定数項を含む）
opt_gauss <- optim(par = rep(1, k+1), fn = LL_gauss, control = list(fnscale = -1), hessian = TRUE) # 尤度最大化
beta_gauss <- opt_gauss$par[1:k] # 偏回帰係数
se_gauss <- diag(sqrt(abs(solve(opt_gauss$hessian[1:k, 1:k])))) # 標準誤差
t_gauss <- beta_gauss / se_gauss # t値
coef_gauss <- cbind(coef = beta_gauss, SE = se_gauss, t = t_gauss)
rss_gauss <- sum((y_gauss - x_gauss %*% beta_gauss)^2) # RSS
tss_gauss <- sum((y_gauss - mean(y_gauss))^2) # TSS
r2_gauss <- 1 - rss_gauss / tss_gauss # 決定係数
aic_gauss <- -2 * opt_gauss$value + 2 * (k+1) # AIC = -2(max LL - #param)
list(coef = coef_gauss, R2 = r2_gauss, RSS = rss_gauss, s2 = opt_gauss$par[k+1], logLik = opt_gauss$value, AIC = aic_gauss)
}計算してみる．
ml_gauss(y_gauss = y_swiss, x_gauss = x_swiss)## Warning in log(2 * pi * param[k + 1]): 計算結果が NaN になりました
## Warning in log(2 * pi * param[k + 1]): 計算結果が NaN になりました## $coef
## coef SE t
## [1,] 94.66923521 7.57484744 12.497841
## [2,] -0.09443892 0.08313466 -1.135975
## [3,] -1.19762789 0.23666545 -5.060426
##
## $R2
## [1] 0.4326027
##
## $RSS
## [1] 4072.752
##
## $s2
## [1] 86.69205
##
## $logLik
## [1] -171.5454
##
## $AIC
## [1] 351.0908summary(lm_swiss3 <- lm(Fertility ~ Agriculture + Examination, swiss)) # 回帰係数の推定量##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.4089 -6.3234 0.0577 6.3134 20.8937
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.6097 7.8271 12.087 1.41e-15 ***
## Agriculture -0.0940 0.0859 -1.094 0.28
## Examination -1.1950 0.2445 -4.887 1.40e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.621 on 44 degrees of freedom
## Multiple R-squared: 0.4326, Adjusted R-squared: 0.4068
## F-statistic: 16.77 on 2 and 44 DF, p-value: 3.85e-06logLik(lm_swiss3) # 対数尤度## 'log Lik.' -171.5453 (df=4)AIC(lm_swiss3)## [1] 351.0906（余談）ところで optim 関数の初期値が先程と異なっています．この程度であれば試行錯誤で何とかなりますが，複雑になるにつれて初期値の設定一つとっても職人芸のようになっていきます…
尤度比検定 (likelihood ratio test) は lmtest::lrtest で計算できます．回帰オブジェクトの対数尤度を返す logLik 関数を用いればパッケージに頼らずとも自分で計算できます．
\[ 2 \ln \left( \frac{L^1}{L^0} \right) = 2 (l^1 - l^0) \sim \chi^2 \]
lm_swiss1 <- lm(Fertility ~ . , swiss)
lm_swiss2 <- lm(Fertility ~ . - Agriculture - Examination, swiss)
# install.packages("lmtest")
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlmtest::lrtest(lm_swiss1, lm_swiss2)## Likelihood ratio test
##
## Model 1: Fertility ~ Agriculture + Examination + Education + Catholic +
## Infant.Mortality
## Model 2: Fertility ~ (Agriculture + Examination + Education + Catholic +
## Infant.Mortality) - Agriculture - Examination
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 7 -156.04
## 2 5 -159.33 -2 6.5969 0.03694 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1llr_swiss <- 2 * (logLik(object = lm_swiss1) - logLik(object = lm_swiss2)) # 対数尤度比
1 - pchisq(q = abs(llr_swiss), df = 2) # p値## 'log Lik.' 0.03694084 (df=7)WLS / 加重最小二乗法
(TBC)
GLM / 一般化線形モデル
古典的正規回帰モデルもGLMの一部として認識できます．
Logit / ロジスティック回帰
See also:
- NaoyaOura 「多変量解析の超基本！知られざるロジスティック回帰分析のアルゴリズム導出」 [Qiita, 2015-06-30]
被説明変数が {0, 1} の場合は二項選択モデル．多項選択の場合は nnet::multinom などを参照されたい．
\[ \rm{logit} (p_i) = \ln \frac{p_i}{1 - p_i} = x_{i}^{'} \beta, \] \[ p_i = F (x_{i}^{'} \beta) = \frac{1}{1 + \exp (- x_{i}^{'} \beta)} \] \[ \ln L = l = \sum_i \left( y_i \ln (F (x_i' \beta)) + (1 - y_i) \ln (1 - F (x_i' \beta)) \right) \]
適当にデータを作って glm 関数によって推定してみましょう．関数の引数として family を指定する点がポイントです．
x2 <- runif(n = 30) # 適当に説明変数を生成 （一様乱数）
x3 <- cbind(1, x2) # 説明変数行列：定数項の 1 を含む
y2 <- round(x2 + runif(n = 30, min = -.5, max = .5)) # 被説明変数
plot(x2, y2) # 描画
summary(glm(formula = y2 ~ x2, family = binomial(link = "logit")))##
## Call:
## glm(formula = y2 ~ x2, family = binomial(link = "logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0667 -0.3936 -0.2465 0.4831 1.9122
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.212 1.434 -2.938 0.00331 **
## x2 8.727 2.834 3.080 0.00207 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.054 on 29 degrees of freedom
## Residual deviance: 21.456 on 28 degrees of freedom
## AIC: 25.456
##
## Number of Fisher Scoring iterations: 5ここまでと同様に，最尤法を用いた推定を行う自作の関数も書いてみましょう．
ml_logit <- function (y_logit, x_logit) {
LL_logit <- function (param) { # 対数尤度関数
plogis_xbeta <- plogis(q = x_logit %*% param, location = 0, scale = 1) # ロジット分布の分布関数
sum(y_logit %*% log(plogis_xbeta) + (1-y_logit) %*% log(1 - plogis_xbeta))
}
k <- ncol(x_logit)
opt_logit <- optim(par = rep(0, k), fn = LL_logit, control = list(fnscale = -1), hessian = TRUE)
se_logit <- diag(sqrt(abs(solve(opt_logit$hessian))))
list(Estimates = cbind(coef = opt_logit$par, SE = se_logit), logLik = opt_logit$value)
}計算してみる．
ml_logit(y_logit = y2, x_logit = x3)## $Estimates
## coef SE
## [1,] -4.212013 1.433764
## [2,] 8.726117 2.833592
##
## $logLik
## [1] -10.72809logLik(glm(formula = y2 ~ x2, family = binomial(link = "logit"))) # 確認## 'log Lik.' -10.72809 (df=2)Probit / プロビット回帰
ロジットとはリンク関数が異なるだけですから，これを標準正規分布 pnorm で設定すればよいだけです．
ml_probit <- function (y_probit, x_probit) {
LL_probit <- function (param) {
pnorm_xbeta <- pnorm(q = x_probit %*% param, mean = 0, sd = 1) # 標準正規分布の分布関数
sum(y_probit %*% log(pnorm_xbeta) + (1-y_probit) %*% log(1 - pnorm_xbeta))
}
optim(par = rep(0, ncol(x_probit)), fn = LL_probit, control = list(fnscale = -1))
}計算してみる．
ml_probit(y_probit = y2, x_probit = x3)## $par
## [1] -2.491804 5.144406
##
## $value
## [1] -10.5941
##
## $counts
## function gradient
## 67 NA
##
## $convergence
## [1] 0
##
## $message
## NULLsummary(glm(formula = y2 ~ x2, family = binomial(link = "probit"))) # 答え合わせ##
## Call:
## glm(formula = y2 ~ x2, family = binomial(link = "probit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0600 -0.3739 -0.1993 0.4845 1.9062
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4926 0.7445 -3.348 0.000814 ***
## x2 5.1458 1.4753 3.488 0.000487 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.054 on 29 degrees of freedom
## Residual deviance: 21.188 on 28 degrees of freedom
## AIC: 25.188
##
## Number of Fisher Scoring iterations: 6optim の―時として不可解な挙動に―信頼をおけないような場合は反復加重最小二乗法（iteratively reweighted least squares, IRLS）で推定することもできます．高野佳祐さんによるコードを RPubs でご確認ください．
Poisson regression / ポアソン回帰
See also: - ポアソン回帰 (奥村先生)
Tobit / トービット
See also: - Rでトービット・モデルによる打ち切りデータの推定
QR / 分位点回帰
Quantile regressionです．quantreg::qr 関数の引数 tau に [0, 1] で分位点を指定すればOK．
# install.packages("quantreg")
library(quantreg)## Loading required package: SparseM##
## Attaching package: 'SparseM'## The following object is masked from 'package:base':
##
## backsolveqr_swiss <- rq(formula = Fertility ~ . , tau = (1:10)/11, data = swiss)
plot(summary(qr_swiss))
Time series / 時系列分析
VAR（ベクトル／多変量自己回帰）モデル．VARは基本統計量の算出や線形回帰モデルの推定よりadvancedな分析のため，Rを立ち上げた時に自動的に読み込まれる標準の関数には含まれていません．そこでVAR分析をカバーする {vars} パッケージをインストールしましょう．
# install.packages("vars") # インストールは最初の1回のみ
library(vars) # 読み込みはRを起動するたびに実行## Warning: package 'vars' was built under R version 3.3.3## Loading required package: MASS## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 3.3.3## Loading required package: sandwich## Loading required package: urca## Warning: package 'urca' was built under R version 3.3.3data(Canada) # データを使用することの宣言：カナダのマクロ経済指標時系列
VARselect(y = Canada[, 1:2], type = "const") # 次数の決定：トレンドなし・定数項あり## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 2 2 3
##
## $criteria
## 1 2 3 4 5
## AIC(n) -2.2547500 -2.79747658 -2.8392941 -2.74087999 -2.66436380
## HQ(n) -2.1802268 -2.67327119 -2.6654066 -2.51731028 -2.39111194
## SC(n) -2.0679339 -2.48611644 -2.4033899 -2.18043172 -1.97937147
## FPE(n) 0.1049091 0.06098885 0.0585333 0.06466987 0.06995344
## 6 7 8 9 10
## AIC(n) -2.58024335 -2.50122558 -2.43048042 -2.36259380 -2.36914450
## HQ(n) -2.25730933 -2.12860941 -2.00818209 -1.89061332 -1.84748186
## SC(n) -1.77070696 -1.56714514 -1.37185592 -1.17942524 -1.06143188
## FPE(n) 0.07631563 0.08292335 0.08947618 0.09641551 0.09660646var_canada <- VAR(y = Canada[, 1:2], p = 2, type = "const") # 推定
summary(var_canada)##
## VAR Estimation Results:
## =========================
## Endogenous variables: e, prod
## Deterministic variables: const
## Sample size: 82
## Log Likelihood: -117.626
## Roots of the characteristic polynomial:
## 0.9962 0.8864 0.6927 0.2621
## Call:
## VAR(y = Canada[, 1:2], p = 2, type = "const")
##
##
## Estimation results for equation e:
## ==================================
## e = e.l1 + prod.l1 + e.l2 + prod.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## e.l1 1.58569 0.08576 18.489 < 2e-16 ***
## prod.l1 0.21176 0.06052 3.499 0.00078 ***
## e.l2 -0.60923 0.08163 -7.463 1.09e-10 ***
## prod.l2 -0.15980 0.06257 -2.554 0.01263 *
## const 1.17425 4.67810 0.251 0.80248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.379 on 77 degrees of freedom
## Multiple R-Squared: 0.9983, Adjusted R-squared: 0.9982
## F-statistic: 1.134e+04 on 4 and 77 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation prod:
## =====================================
## prod = e.l1 + prod.l1 + e.l2 + prod.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## e.l1 0.2494 0.1566 1.593 0.11530
## prod.l1 1.2516 0.1105 11.327 < 2e-16 ***
## e.l2 -0.2200 0.1490 -1.476 0.14409
## prod.l2 -0.3208 0.1142 -2.808 0.00631 **
## const 0.4279 8.5414 0.050 0.96017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.692 on 77 degrees of freedom
## Multiple R-Squared: 0.9747, Adjusted R-squared: 0.9734
## F-statistic: 742.5 on 4 and 77 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## e prod
## e 0.14364 0.01694
## prod 0.01694 0.47884
##
## Correlation matrix of residuals:
## e prod
## e 1.00000 0.06458
## prod 0.06458 1.00000plot(var_canada)
serial.test(x = var_canada) # 残差に対する系列相関の検定##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var_canada
## Chi-squared = 48.188, df = 56, p-value = 0.7617normality.test(x = var_canada) # 残差に対する正規性の検定## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object var_canada
## Chi-squared = 3.7372, df = 4, p-value = 0.4427
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object var_canada
## Chi-squared = 0.73129, df = 2, p-value = 0.6937
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object var_canada
## Chi-squared = 3.0059, df = 2, p-value = 0.2225arch.test(x = var_canada) # 残差に対するARCH-LM検定##
## ARCH (multivariate)
##
## data: Residuals of VAR object var_canada
## Chi-squared = 59.767, df = 45, p-value = 0.06923Granger因果性検定．
causality(x = var_canada, cause = "e") # e: employment## $Granger
##
## Granger causality H0: e do not Granger-cause prod
##
## data: VAR object var_canada
## F-Test = 1.7025, df1 = 2, df2 = 154, p-value = 0.1856
##
##
## $Instant
##
## H0: No instantaneous causality between: e and prod
##
## data: VAR object var_canada
## Chi-squared = 0.3406, df = 1, p-value = 0.5595# causality(x = var_canada, cause = "prod") # prod: labour productivity
# causality(x = var_canada, cause = "rw") # rw: real wage
# causality(x = var_canada, cause = "U") # U: unemployment rateインパルス応答関数．
irf_canada <- irf(x = var_canada, n.ahead = 20)
plot(irf_canada)
Panel analysis
- 太郎丸博: Rでパネルデータ分析
通常のパネル分析．米国の Gross State Product (gsp) の分析を例にして．
library(plm)## Loading required package: Formuladata(Produc, package = "plm") # [states] x [1970-1986]
fm_produc <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
# Pooling分析は，全データをcross-sectionデータのように分析した場合と同じ
p_produc_pooling <- plm(formula = fm_produc, data = Produc, model = "pooling") # pooled OLS
produc_pooling <- lm(formula = fm_produc, data = Produc)
# 固定効果モデル：個体FE
p_produc_FE <- plm(formula = fm_produc, data = Produc, model = "within")
summary(p_produc_FE)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = fm_produc, data = Produc, model = "within")
##
## Balanced Panel: n=48, T=17, N=816
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.12000 -0.02370 -0.00204 0.01810 0.17500
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## log(pcap) -0.02614965 0.02900158 -0.9017 0.3675
## log(pc) 0.29200693 0.02511967 11.6246 < 2.2e-16 ***
## log(emp) 0.76815947 0.03009174 25.5273 < 2.2e-16 ***
## unemp -0.00529774 0.00098873 -5.3582 1.114e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 18.941
## Residual Sum of Squares: 1.1112
## R-Squared: 0.94134
## Adj. R-Squared: 0.93742
## F-statistic: 3064.81 on 4 and 764 DF, p-value: < 2.22e-16produc_FE <- lm(formula = update.formula(old = fm_produc, new = . ~ . + state), data = Produc) # cf. LSDV
head(summary(plm::fixef(p_produc_FE, type = "dmean")), 2) # dmean: deviation from the overall mean## Estimate Std. Error t-value Pr(>|t|)
## ALABAMA -0.15128205 0.1760039 -0.85953818 0.3903133
## ARIZONA 0.01518903 0.1751885 0.08670108 0.9309318# 時間FE
p_produc_FEtime <- plm(formula = fm_produc, data = Produc, effect = "time", model = "within")
head(summary(plm::fixef(p_produc_FEtime, type = "dmean")), 2)## Estimate Std. Error t-value Pr(>|t|)
## 1970 -0.008568690 0.05814126 -0.14737710 0.8828717
## 1971 -0.002577321 0.05850905 -0.04404996 0.9648756# 個体FEと時間FEの両方
p_produc_FEtwo <- plm(formula = fm_produc, data = Produc, effect = "twoways", model = "within")
# summary(plm::fixef(p_produc_FEtwo, type = "dmean")) # for individual (default)
# summary(plm::fixef(p_produc_FEtwo, effect = "time", type = "dmean"))
# ランダム効果モデル
p_produc_RE <- plm(formula = fm_produc, data = Produc, model = "random") # Random effect model
phtest(x = p_produc_FE, x2 = p_produc_RE) # Hausman Test (FE vs. RE)##
## Hausman Test
##
## data: fm_produc
## chisq = 9.5254, df = 4, p-value = 0.04923
## alternative hypothesis: one model is inconsistent留意すべき点ですが，plm::plm 関数では data （data.frame型）の1列目が individual，2列目が time を示す変数とみなされます．以下の例で確認してみてください．
Produc2 <- cbind(year = Produc[, 2], state = Produc[, 1], Produc[, 3:11])
p_produc_FE2 <- plm(formula = fm_produc, data = Produc2)
summary(p_produc_FE2) # summary(p_produc_FEtime) と同じ## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = fm_produc, data = Produc2)
##
## Balanced Panel: n=17, T=48, N=816
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.22400 -0.05830 -0.00135 0.04900 0.35800
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## log(pcap) 0.1647800 0.0174912 9.4207 < 2.2e-16 ***
## log(pc) 0.3035960 0.0104427 29.0727 < 2.2e-16 ***
## log(emp) 0.5888107 0.0137757 42.7428 < 2.2e-16 ***
## unemp -0.0060575 0.0017702 -3.4220 0.0006534 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 834.46
## Residual Sum of Squares: 6.0629
## R-Squared: 0.99273
## Adj. R-Squared: 0.99255
## F-statistic: 27156 on 4 and 795 DF, p-value: < 2.22e-16Spatial regression
Spatial dependence
SWM / 空間重み行列
まず空間重み行列（SWM）を計算します．使用するデータ columbus はオハイオ州コロンバスの犯罪に関連するデータで，パッケージ {spdep} で見本データとして用意されています．
\[ w_{ij} = 1 / d_{ij}^{2} \]
# install.packages("spdep")
library(spdep)## Loading required package: sp## Loading required package: Matrixdata(columbus, package = "spdep") # Columbus OH spatial analysis data set
income7 <- cut(x = columbus$INC, breaks = 7) # INC (income) の値を幅の等しい幾つかの区間に分割
plot(columbus$X, columbus$Y, col = rainbow(n = 8)[income7]) # 区間ごとに色を変えて描画
legend(x = "topleft", legend = levels(income7), col = rainbow(8)[1:7], pch = 1) # 凡例
dist_col <- as.matrix(dist(columbus[, c("X", "Y")])) # 距離行列
dist_col[dist_col > 5] <- Inf # cutt-off (距離 5 で cut)
w_col <- 1 / dist_col^2 # 重力モデルのアナロジーより alpha = 2
diag(w_col) <- 0 # 対角要素は 0
w_std_col <- w_col / apply(X = w_col, MARGIN = 1, FUN = sum) # 行基準化
listw_std_col <- mat2listw(x = w_std_col) # matrix 型 → list 型重み行列が計算できればモランのIが計算できます．
\[ I = \frac{n}{\sum_i \sum_j w_{ij}} \frac{e' W e}{e' e} \]
fm_col <- CRIME ~ INC + HOVAL
lm_col <- lm(formula = fm_col, data = columbus) # OLS
resid_col <- residuals(lm_col) # 回帰残差
resid_col %*% w_std_col %*% resid_col / crossprod(resid_col) # W は行基準化されているため n = ΣΣw## [,1]
## [1,] 0.2733901moran.test(x = resid_col, listw = listw_std_col) # Moran's I test: 答え合わせ##
## Moran I test under randomisation
##
## data: resid_col
## weights: listw_std_col
##
## Moran I statistic standard deviate = 3.4369, p-value = 0.0002942
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.27339015 -0.02083333 0.00732875# 空間的自己相関の診断
lm.LMtests(model = lm_col, listw = listw_std_col, test = "all")## Warning in lm.LMtests(model = lm_col, listw = listw_std_col, test = "all"):
## Spatial weights matrix not row standardized##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = fm_col, data = columbus)
## weights: listw_std_col
##
## LMerr = 8.9376, df = 1, p-value = 0.002794
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = fm_col, data = columbus)
## weights: listw_std_col
##
## LMlag = 15.11, df = 1, p-value = 0.0001014
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = fm_col, data = columbus)
## weights: listw_std_col
##
## RLMerr = 0.10569, df = 1, p-value = 0.7451
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = fm_col, data = columbus)
## weights: listw_std_col
##
## RLMlag = 6.2778, df = 1, p-value = 0.01223
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = fm_col, data = columbus)
## weights: listw_std_col
##
## SARMA = 15.215, df = 2, p-value = 0.0004966{spdep} パッケージで定義されている関数を使って重み行列を計算することもできます．
coord_col <- coordinates(columbus[, c("X", "Y")]) # as.matrixでもOK
dnb_col <- dnearneigh(x = coord_col, d1 = 0, d2 = 5) # [0, 5] の距離にある個票のペア
dist_list_col <- nbdists(nb = dnb_col, coords = coord_col) # ペア間の距離
nb_w_list_col <- lapply(X = dist_list_col, FUN = function (x) 1 / x^2)
w_list_col <- nb2listw(neighbours = dnb_col, glist = nb_w_list_col, style = "W") # list型のSWM
w_mat_col <- nb2mat(neighbours = dnb_col, glist = nb_w_list_col, style = "W") # matrix型のSWM
moran.test(x = resid_col, listw = w_list_col) # Moran's I##
## Moran I test under randomisation
##
## data: resid_col
## weights: w_list_col
##
## Moran I statistic standard deviate = 3.4369, p-value = 0.0002942
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.27339015 -0.02083333 0.00732875Spatial econometrics / 空間計量経済モデル
空間誤差モデル (SEM)，空間階差モデル (SLM)，空間ダービンモデル (SDM)，空間的自己相関を考慮した条件付き自己回帰 (CAR) モデルも以下のように計算できます（全て最尤法で推定）．
SEM: \[ y = X \beta + u, u = \lambda W u + \varepsilon, \varepsilon \sim N (0, \sigma^2 I) \] SLM: \[ y = \rho W y + X \beta + \varepsilon, \varepsilon \sim N (0, \sigma^2 I) \] SDM: \[ y = \rho W y + X \beta + W X \beta \gamma + \varepsilon, \varepsilon \sim N (0, \sigma^2 I) \]
sem_col <- errorsarlm(formula = fm_col, data = columbus, listw = w_list_col) # SEM
summary(sem_col)##
## Call:errorsarlm(formula = fm_col, data = columbus, listw = w_list_col)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.60925 -6.61352 -0.10169 5.90601 22.96155
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 54.663185 5.882899 9.2919 < 2.2e-16
## INC -0.870672 0.312933 -2.7823 0.005398
## HOVAL -0.256044 0.084609 -3.0262 0.002476
##
## Lambda: 0.65634, LR test value: 12.315, p-value: 0.00044924
## Asymptotic standard error: 0.11862
## z-value: 5.5333, p-value: 3.1421e-08
## Wald statistic: 30.618, p-value: 3.1421e-08
##
## Log likelihood: -181.2195 for error model
## ML residual variance (sigma squared): 85.166, (sigma: 9.2285)
## Number of observations: 49
## Number of parameters estimated: 5
## AIC: 372.44, (AIC for lm: 382.75)moran.test(x = residuals(sem_col), listw = w_list_col)##
## Moran I test under randomisation
##
## data: residuals(sem_col)
## weights: w_list_col
##
## Moran I statistic standard deviate = 0.92753, p-value = 0.1768
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.057232866 -0.020833333 0.007083823slm_col <- lagsarlm(formula = fm_col, data = columbus, listw = w_list_col) # SLM
sdm_col <- lagsarlm(formula = fm_col, data = columbus, listw = w_list_col, type = "mixed") # SDM
car_col <- spautolm(formula = fm_col, data = columbus, listw = w_list_col, family = "CAR") # CAR## Warning in spautolm(formula = fm_col, data = columbus, listw =
## w_list_col, : Non-symmetric spatial weights in CAR model## Warning in sqrt(fdHess[1, 1]): 計算結果が NaN になりましたSEMを errorsarlm 関数に頼らないで回すには for 文による繰り返し計算のコードを自作する必要があります．研究室の後輩である高野佳祐さんによるコードを RPubs に転載しましたので，こちらをご覧ください．
Spatial panel analysis / 空間パネル分析
library(splm)## Warning: package 'splm' was built under R version 3.3.3data(Produc, package = "plm")
data(usaww) # SWM for US states
fm_produc <- log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
# FE panel with spatial errors
sp_produc_FE <- spml(formula = fm_produc, data = Produc, listw = mat2listw(usaww))
summary(sp_produc_FE)## Spatial panel fixed effects error model
##
##
## Call:
## spml(formula = fm_produc, data = Produc, listw = mat2listw(usaww))
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.1250 -0.0238 -0.0035 0.0171 0.1880
##
## Spatial error parameter:
## Estimate Std. Error t-value Pr(>|t|)
## rho 0.557401 0.033075 16.853 < 2.2e-16 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## log(pcap) 0.0051438 0.0250109 0.2057 0.83705
## log(pc) 0.2053026 0.0231427 8.8712 < 2e-16 ***
## log(emp) 0.7822540 0.0278057 28.1328 < 2e-16 ***
## unemp -0.0022317 0.0010709 -2.0839 0.03717 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1effects.splm(sp_produc_FE)##
## Intercept:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 2.84695 0.14324 19.876 < 2.2e-16 ***
##
##
## Spatial fixed effects:
## Estimate Std. Error t-value Pr(>|t|)
## 1 -0.1393449 0.1442142 -0.9662 0.33393
## 2 0.0042234 0.1435461 0.0294 0.97653
## 3 -0.0993829 0.1369773 -0.7255 0.46812
## 4 0.2059001 0.1648750 1.2488 0.21173
## 5 0.0422873 0.1417953 0.2982 0.76553
## 6 0.0912363 0.1410148 0.6470 0.51763
## 7 -0.0168266 0.1325054 -0.1270 0.89895
## 8 0.0084710 0.1501192 0.0564 0.95500
## 9 -0.0791192 0.1452860 -0.5446 0.58605
## 10 -0.0556094 0.1308325 -0.4250 0.67081
## 11 0.0854043 0.1557281 0.5484 0.58340
## 12 -0.0570390 0.1455946 -0.3918 0.69523
## 13 -0.0121360 0.1468313 -0.0827 0.93413
## 14 0.0088455 0.1454997 0.0608 0.95152
## 15 0.0459446 0.1466921 0.3132 0.75413
## 16 0.2339888 0.1535306 1.5241 0.12750
## 17 -0.1393260 0.1269724 -1.0973 0.27251
## 18 0.0284226 0.1471209 0.1932 0.84681
## 19 0.0017597 0.1437518 0.0122 0.99023
## 20 0.0915863 0.1533967 0.5971 0.55047
## 21 -0.0265202 0.1498833 -0.1769 0.85956
## 22 -0.1163081 0.1408047 -0.8260 0.40879
## 23 -0.0232363 0.1446846 -0.1606 0.87241
## 24 -0.0194968 0.1415010 -0.1378 0.89041
## 25 -0.0556076 0.1470265 -0.3782 0.70527
## 26 -0.0384417 0.1328683 -0.2893 0.77234
## 27 -0.0984937 0.1247750 -0.7894 0.42990
## 28 0.1074551 0.1456143 0.7379 0.46055
## 29 0.0987881 0.1393024 0.7092 0.47822
## 30 0.1565340 0.1648349 0.9496 0.34229
## 31 -0.0944648 0.1419280 -0.6656 0.50568
## 32 -0.0178559 0.1407789 -0.1268 0.89907
## 33 0.0127804 0.1538395 0.0831 0.93379
## 34 0.0998038 0.1437757 0.6942 0.48758
## 35 -0.0189651 0.1411083 -0.1344 0.89309
## 36 -0.0229365 0.1546517 -0.1483 0.88210
## 37 -0.0626022 0.1230998 -0.5085 0.61107
## 38 -0.2308038 0.1359967 -1.6971 0.08967 .
## 39 -0.1039776 0.1384520 -0.7510 0.45265
## 40 -0.1178466 0.1471952 -0.8006 0.42335
## 41 0.1504499 0.1608640 0.9353 0.34965
## 42 -0.0647586 0.1366046 -0.4741 0.63546
## 43 -0.1132443 0.1254047 -0.9030 0.36651
## 44 0.0311400 0.1458533 0.2135 0.83094
## 45 0.0864727 0.1541201 0.5611 0.57475
## 46 -0.0594311 0.1405449 -0.4229 0.67240
## 47 -0.0215054 0.1467946 -0.1465 0.88353
## 48 0.3137863 0.1466038 2.1404 0.03232 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# RE panel with spatial lag
sp_produc_RE <- spml(formula = fm_produc, data = Produc, listw = mat2listw(usaww), model = "random", spatial.error = "none", lag = TRUE)
summary(sp_produc_RE)## ML panel with spatial lag, random effects
##
## Call:
## spreml(formula = formula, data = data, index = index, w = listw2mat(listw),
## w2 = listw2mat(listw2), lag = lag, errors = errors, cl = cl)
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.38 1.57 1.70 1.70 1.80 2.13
##
## Error variance parameters:
## Estimate Std. Error t-value Pr(>|t|)
## phi 21.3175 8.3017 2.5678 0.01023 *
##
## Spatial autoregressive coefficient:
## Estimate Std. Error t-value Pr(>|t|)
## lambda 0.161615 0.029099 5.554 2.793e-08 ***
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 1.65814995 0.15071855 11.0016 < 2.2e-16 ***
## log(pcap) 0.01294505 0.02493997 0.5190 0.6037
## log(pc) 0.22555376 0.02163422 10.4258 < 2.2e-16 ***
## log(emp) 0.67081075 0.02642113 25.3892 < 2.2e-16 ***
## unemp -0.00579716 0.00089175 -6.5009 7.984e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Spatially varying coefficient model: GWR
地理的加重回帰 (GWR) 分析．偏回帰係数パラメータに（地理的な意味で）局所的な多重共線性が発生するという問題から実用性については謎に包まれているようです．
\[ W_i^{1/2} y = W_i^{1/2} X \beta_i + W_i^{1/2} \varepsilon \] \[ w_{ij} = \exp \left( - \frac{d_{ij}^{2}}{2 \theta^2} \right) \] \[ W_i = \left( {\begin{array} w_{i1} & \cdots & 0 \\ \vdots & w_{ij} & \vdots \\ 0 & \cdots & w_{in} \\ \end{array}} \right) \] \[ \hat{\theta} = \arg \min \sum_i (y_i - \hat{y}_{\ne i} (\theta))^2 \] \[ \hat{\beta_i} = (X' W_i X)^{-1} X' W_i y \]
# install.packages("spgwr")
library(spgwr)## Warning: package 'spgwr' was built under R version 3.3.3## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)data(columbus, package = "spgwr")
fm_col <- crime ~ income + housing
coords_col <- cbind(columbus$x, columbus$y) # 座標
bw_col <- gwr.sel(formula = fm_col, data = columbus, coords = coords_col) # バンド幅を CV で最適化## Bandwidth: 12.65221 CV score: 7432.171
## Bandwidth: 20.45127 CV score: 7462.669
## Bandwidth: 7.832129 CV score: 7323.494
## Bandwidth: 4.853154 CV score: 7307.484
## Bandwidth: 5.122256 CV score: 7322.64
## Bandwidth: 3.012046 CV score: 6461.665
## Bandwidth: 1.874178 CV score: 6473.183
## Bandwidth: 2.475225 CV score: 6109.765
## Bandwidth: 2.447682 CV score: 6098.241
## Bandwidth: 2.228623 CV score: 6063.971
## Bandwidth: 2.264529 CV score: 6060.644
## Bandwidth: 2.280636 CV score: 6060.52
## Bandwidth: 2.274941 CV score: 6060.473
## Bandwidth: 2.275073 CV score: 6060.473
## Bandwidth: 2.275032 CV score: 6060.473
## Bandwidth: 2.274991 CV score: 6060.473
## Bandwidth: 2.275032 CV score: 6060.473gwr_col <- gwr(formula = fm_col, data = columbus, coords = coords_col, bandwidth = bw_col, hatmatrix = TRUE)
gwr_col## Call:
## gwr(formula = fm_col, data = columbus, coords = coords_col, bandwidth = bw_col,
## hatmatrix = TRUE)
## Kernel function: gwr.Gauss
## Fixed bandwidth: 2.275032
## Summary of GWR coefficient estimates at data points:
## Min. 1st Qu. Median 3rd Qu. Max. Global
## X.Intercept. 23.23000 54.13000 63.90000 68.76000 80.90000 68.6189
## income -3.13100 -1.91300 -0.98440 -0.36860 1.29100 -1.5973
## housing -1.05300 -0.37670 -0.09739 0.03006 0.79460 -0.2739
## Number of data points: 49
## Effective number of parameters (residual: 2traceS - traceS'S): 29.61664
## Effective degrees of freedom (residual: 2traceS - traceS'S): 19.38336
## Sigma (residual: 2traceS - traceS'S): 8.027396
## Effective number of parameters (model: traceS): 23.92826
## Effective degrees of freedom (model: traceS): 25.07174
## Sigma (model: traceS): 7.058251
## Sigma (ML): 5.048836
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 403.6193
## AIC (GWR p. 96, eq. 4.22): 321.6617
## Residual sum of squares: 1249.046
## Quasi-global R2: 0.9070521sp::spplot(obj = gwr_col$SDF, zcol = "income") # "income" の回帰係数
# sp::spplot(obj = gwr_col$SDF, zcol = "localR2") # local R2簡単のためにバンド幅が所与の元での地点ごとの回帰係数を推定してみましょう．
gwr2 <- function (y_gwr, x_gwr, coords_gwr, bw_gwr) {
# y_gwr 被説明変数ベクトル
# x_gwr 説明変数行列
# coords_gwr 座標行列（n x 2）
# bw_gwr バンド幅（この関数では外生で与える）
beta_mat <- NULL
dist_gwr <- as.matrix(dist(coords_gwr))
for (i in 1:length(y_gwr)) {
W_i <- diag(exp(- .5 * (dist_gwr[i, ] / bw_gwr)^2)) # Gaussian kernel: 「gwr.Gauss」に対応
beta_i <- as.numeric(solve(t(x_gwr) %*% W_i %*% x_gwr) %*% t(x_gwr) %*% W_i %*% y_gwr)
e_i <- y_gwr[i] - t(x_gwr[i, ]) %*% beta_i # 残差
beta_mat <- rbind(beta_mat, c(beta_i, e_i))
}
beta_mat
}計算してみる．
x_col <- cbind(1, columbus$income, columbus$housing) # X
gwr_col2 <- gwr2(y_gwr = columbus$crime, x_gwr = x_col, coords_gwr = coords_col, bw_gwr = 2.275)
plot(gwr_col$SDF@data$income, gwr_col2[,2]) # spgwr::gwr の推定結果と比較
abline(a = 0, b = 1) # 45度線
Kriging / クリギング
gstat パッケージを使えば簡単にkrigingが実行できます．ここでは regression kriging （universal krigingと表記されることが多いかも）を例示します．
# install.packages("gstat")
library(gstat)
data(meuse, package = "sp") # 欧州・マース川沿いにおける重金属データ
# str(meuse)
# plot(x = meuse$dist, y = log(meuse$zinc))
obj_m <- log(zinc) ~ x + y + dist # Zn: 亜鉛
loc_m <- ~ x + y # 座標を表す変数（データフレーム中の列名）の指定
vario_m <- variogram(object = obj_m, locations = loc_m, data = meuse, cressie = T)
# plot(vario_m) # バリオグラム
# plot(variogram(obj_m, loc_m, meuse, cressie = T, cloud = T)) # バリオグラム雲
vgm_m <- vgm(psill = 0.2, model = "Sph", range = 800, nugget = 0.1) # 球形のバリオグラム・モデル
fitvario_m <- fit.variogram(object = vario_m, model = vgm_m) # 当てはまりの良いパラメータに調整
plot(vario_m, fitvario_m)
バリオグラム・モデルを特定化しパラメータを決定したため，いよいよ内挿を行います．
gstat_m <- gstat(id = "zn", formula = obj_m, locations = loc_m, data = meuse, model = fitvario_m)
data(meuse.grid, package = "sp") # 説明変数のみ
pred_m <- predict(object = gstat_m, newdata = meuse.grid) # 内挿## [using universal kriging]head(pred_m)## x y zn.pred zn.var
## 1 181180 333740 6.803302 0.1869643
## 2 181140 333700 6.839184 0.1573393
## 3 181180 333700 6.744413 0.1655273
## 4 181220 333700 6.601696 0.1743819
## 5 181100 333660 6.884581 0.1270611
## 6 181140 333660 6.783586 0.1355028pred_m7 <- cut(x = pred_m$zn.pred, breaks = 7) # 等幅区間に分割
plot(x = pred_m$x, y = pred_m$y, col = terrain.colors(n = 7)[pred_m7], pch = 20) # 描画
legend(x = "topleft", legend = levels(pred_m7), col = terrain.colors(7), pch = 20) # 凡例
ヘルプ（?gstat）に従って lattice::levelplot を用いれば上記より多少マシなコロプレス図を得ることができます．
Hypothesis testing / 統計的仮説検定
おまけ．
統計的仮説検定 (statistical hypothesis testing; null hypothesis significance testing = NHST) は，近年（特に心理学など一部の分野では）敬遠される傾向もあるようです．HARKing (hypothesizing after the results are Known)やp-hackingにも留意する必要があります．
- 竹中明夫: 「統計学的な有意性検定の意味のなさ」
- Chris Woolston: Psychology journal bans P values, Nature, 519 (7541), Feb 2015.
- スネハ・クルカルニ: 「T&F社のBASP誌がP値の使用禁止を発表」 エディテージ, 2015年10月
とは言っても実証エコノメの世界から仮説検定が一掃されることは当分なさそうですから，各種検定についても一応触れざるを得ません．
Normality test / 正規性の検定
帰無仮説「標本は正規母集団からサンプリングされた」に対して p 値が十分に小さければ（たとえば 有意水準 0.05 より小さければ），H0は棄却（≒正規分布に従っていないだろうと判断）されます．
rand1 <- rnorm(n = 30) # 適当な正規乱数
shapiro.test(x = rand1) # Shapiro-Wilk normality test##
## Shapiro-Wilk normality test
##
## data: rand1
## W = 0.94815, p-value = 0.1508ks.test(x = rand1, y = "pnorm") # Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: rand1
## D = 0.18387, p-value = 0.2322
## alternative hypothesis: two-sidedtseries::jarque.bera.test(x = rand1) # Jarque Bera test##
## Jarque Bera Test
##
## data: rand1
## X-squared = 7.5142, df = 2, p-value = 0.02335Two sample t-test / 二標本の平均値の差の検定
以下のように関数 t.test を実行すると Welch Two Sample t-test と表示され，デフォルトの設定（＝任意指定の引数をいじらず規定値のままにしている場合）では自動的にウェルチのt検定が行われます（等分散であることが既知ならば通常のt検定を実行するよう var.equal 引数を修正すればよい）．
rand1 <- rnorm(n = 30)
rand2 <- runif(n = 10) # 適当な一様乱数
t.test(x = rand1, y = rand2)##
## Welch Two Sample t-test
##
## data: rand1 and rand2
## t = -3.1964, df = 37.51, p-value = 0.002823
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.1137302 -0.2497876
## sample estimates:
## mean of x mean of y
## -0.2114010 0.4703579等分散であるか否かは var.test などを用いれば容易に検定できますが，その結果次第でウェルチの方法による自由度の補正を行うかどうかを決定するのではなく，常時ウェルチの方法を用いてよいでしょう（多重検定の問題が生じるとの指摘もある）．一応，通常のt検定を行う関数を自作してみましょう．
t.test2 <- function (var1, var2) {
n1 <- length(var1)
n2 <- length(var2)
dof <- n1 + n2 - 2 # degree of freedom
var12 <- ((n1-1)*var(var1) + (n2-1)*var(var2)) / dof
t_value <- (mean(var1) - mean(var2)) / sqrt(var12 * (1/n1 + 1/n2))
p_value <- 2*(1 - pt(q = abs(t_value), df = dof))
list(t = t_value, df = dof, p = p_value)
}
t.test2(var1 = rand1, var2 = rand2)## $t
## [1] -2.00662
##
## $df
## [1] 38
##
## $p
## [1] 0.05194432t.test(x = rand1, y = rand2, var.equal = T) # 答え合わせ##
## Two Sample t-test
##
## data: rand1 and rand2
## t = -2.0066, df = 38, p-value = 0.05194
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.369556855 0.006038968
## sample estimates:
## mean of x mean of y
## -0.2114010 0.4703579今更ですが，上記のt検定は二標本が不等分散であろうともそれぞれが正規母集団からのサンプリングであることを仮定していました．この正規性の仮定が満たされない場合は distribution free な non-parametric 検定を行うのが定石です．
wilcox.test(x = rand1, y = rand2) # Wilcoxon rank sum test##
## Wilcoxon rank sum test
##
## data: rand1 and rand2
## W = 79, p-value = 0.02596
## alternative hypothesis: true location shift is not equal to 0ここまでは二標本が独立であることを前提としていましたが，対応がある場合（たとえば，各個体に対して treatment 前後で何らかの測定が行われた場合）はどうでしょうか？
rand3 <- rnorm(n = 10)
rand4 <- runif(n = 10)
t.test(x = rand3, y = rand4, paired = TRUE) # Paired t-test##
## Paired t-test
##
## data: rand3 and rand4
## t = -0.32055, df = 9, p-value = 0.7559
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.8134120 0.6114978
## sample estimates:
## mean of the differences
## -0.1009571wilcox.test(x = rand3, y = rand4, paired = TRUE) # Wilcoxon signed rank test##
## Wilcoxon signed rank test
##
## data: rand3 and rand4
## V = 20, p-value = 0.4922
## alternative hypothesis: true location shift is not equal to 0ウィルコクソンの符号順位検定 (Wilcoxon signed-rank test) ≠ ウィルコクソンの順位和検定 (Wilcoxon rank sum test) ＝ マン・ホイットニーのU検定 (Mann-Whitney U test) です．
（以上）