専門パッケージに頼らないGLM・空間計量経済モデル推定のコーディング
Keisuke Takano (and Sho Kuroda)
May 3, 2017
高野佳佑さん（筑波大学大学院）によって書かれたコードです（黒田は若干調整している程度）．詳細は高野さんにお問い合わせください．spdepなどに頼らず空間計量のモデルの推定を試みているとのことです．デフォルトで含まれているパッケージ群の諸関数には頼っています．
GLM
被説明変数 y が二値変数やカウントデータであった場合．
Logit and probit
反復加重最小二乗法（iteratively reweighted least squares, IRLS）で推定してみます．ロジットもプロビットもプログラムとしてはほとんど同じですので，引数 typeLogitによってロジットかプロビットかを指定し，それに応じて関数の一部が変わるようになっています．
\[ \eta = X \beta \] \[ \mu = \frac{1}{1 + \exp(- \eta)} = \frac{1}{1 + \exp(- X \beta)}, z = \eta + \frac{y - \mu}{\mu (1 - \mu)}, \rm{diag}(W) = w = \mu (1 - \mu) \] \[ \phi = \frac{\exp(-\eta^2 / 2)}{\sqrt{2 \pi}}, \mu = \Phi(\eta) = \Phi(X \beta), z = \frac{\eta + (y - \mu)}{\phi}, \rm{diag}(W) = w = \frac{\phi^2}{\mu (1 - \mu)} \] \[ \beta = \frac{X' W z}{X' W X} \]
irls_logit <- function(y_irls, X_irls, typeLogit = TRUE, maxit = 100, tol = 1e-08) {
# y_irls 被説明変数ベクトル
# x_irls 説明変数行列
# typeLogit TRUEはlogit，FALSEはprobitを分布として仮定
# maxit （収束しない場合の）最大iteration回数
# tol 収束判定に利用するtolerance
y <- as.vector(y_irls)
X <- as.matrix(X_irls)
reg.X <- cbind(1, X_irls)
beta_irls <- rep(0, ncol(reg.X)) # 回帰係数の初期値
for (i in 1:maxit) {
eta <- reg.X %*% beta_irls # 線形予測子
if (typeLogit) { # logit
mu <- 1 / (1 + exp(-eta)) # 平均項
w <- mu * (1-mu)
z <- eta + (y - mu) / (mu * (1-mu))
} else { # probit
phi <- exp((-(eta)^2)/2) / sqrt(2*pi)
mu <- pnorm(q = eta, mean = 0, sd = 1) # 平均項
w <- (phi^2) / (mu*(1-mu))
z <- eta + (y - mu) / phi
}
W <- diag(nrow(reg.X))
diag(W) <- w
old.beta <- beta_irls
beta_irls <- solve(crossprod(reg.X, W %*% reg.X)) %*% crossprod(reg.X, W %*% z)
if(sqrt(crossprod(beta_irls - old.beta)) < tol) break # 収束したら終了
}
beta1 <- as.vector(beta_irls) # 回帰係数
var.beta <- solve(crossprod(reg.X, W %*% reg.X))
stderr <- sqrt(diag(var.beta)) # 回帰係数の標準誤差
z.value <- beta1 / stderr # z値
coef <- cbind(coef = beta1, SE = stderr, z = z.value)
list(coefficient = coef, iteration = i)
}計算してみる．
# 外部サイトのサンプルデータ
bin <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
y_bin <- as.vector(bin[, 1]) # n by 1, {1, 0}
X_bin <- as.matrix(bin[, 2:3]) # n by (k-1)
# IRLS
irls_logit(y_irls = y_bin, X_irls = X_bin) # logit## $coefficient
## coef SE z
## -4.949378063 1.075093072 -4.603674
## gre 0.002690684 0.001057491 2.544403
## gpa 0.754686856 0.319585633 2.361454
##
## $iteration
## [1] 5irls_logit(y_irls = y_bin, X_irls = X_bin, typeLogit = F) # probit## $coefficient
## coef SE z
## -3.003536085 0.6339726699 -4.737643
## gre 0.001642537 0.0006333341 2.593476
## gpa 0.454574944 0.1915197143 2.373515
##
## $iteration
## [1] 7# 答え合わせ
summary(glm(admit ~ gre + gpa, binomial(link = "logit"), bin))##
## Call:
## glm(formula = admit ~ gre + gpa, family = binomial(link = "logit"),
## data = bin)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2730 -0.8988 -0.7206 1.3013 2.0620
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.949378 1.075093 -4.604 4.15e-06 ***
## gre 0.002691 0.001057 2.544 0.0109 *
## gpa 0.754687 0.319586 2.361 0.0182 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 480.34 on 397 degrees of freedom
## AIC: 486.34
##
## Number of Fisher Scoring iterations: 4summary(glm(admit ~ gre + gpa, binomial(link = "probit"), bin))##
## Call:
## glm(formula = admit ~ gre + gpa, family = binomial(link = "probit"),
## data = bin)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2659 -0.9025 -0.7212 1.3018 2.0956
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0035357 0.6339727 -4.738 2.16e-06 ***
## gre 0.0016425 0.0006333 2.593 0.0095 **
## gpa 0.4545748 0.1915198 2.374 0.0176 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 480.19 on 397 degrees of freedom
## AIC: 486.19
##
## Number of Fisher Scoring iterations: 4Regression models for count data: Poisson regression
ポアソン回帰モデルの推定に用いるのは反復加重最小二乗法（iteratively reweighted least squares, IRLS. See ドブソン, 2008: amazon.co.jp）です．
irls_poisson <- function(y_poi, X_poi, maxit = 25, tol = 1e-08, offset = 0) {
# maxit （収束しない場合の）iteration回数
# tol 収束判定に用いるtolerance
y_poi <- as.vector(y_poi) # 被説明変数
X_poi <- as.matrix(X_poi) # 説明変数
reg.X <- cbind(1, X_poi) # +定数項
beta_poi <- rep(0, ncol(reg.X)) # 回帰係数の初期値を設定
for(i in 1:maxit){
eta <- reg.X %*% beta_poi + offset # 線形予測子
mu <- exp(eta) # 平均項
W <- diag(nrow(reg.X)) # 単位行列
diag(W) <- mu # 重み行列（対角項に平均 mu を代入）
z <- eta + (y_poi - mu)*(1/mu) - offset # 調整被説明変数(adjusted dependent variable)
old.beta <- beta_poi
beta_poi <- solve(crossprod(reg.X, W %*% reg.X)) %*% crossprod(reg.X, W %*% z) # 回帰係数を更新
if(sqrt(crossprod(beta_poi - old.beta)) < tol) break # 収束したら終了
}
beta1 <- as.vector(beta_poi) # 回帰係数
var.beta <- solve(crossprod(reg.X, W %*% reg.X)) # 回帰係数の分散
stderr <- sqrt(diag(var.beta)) # 標準誤差
z.value <- beta1 / stderr # z値
# 結果の出力
coef <- cbind(coef = beta1, SE = stderr, z = z.value)
list(coefficient = coef, iteration = i)
}計算してみる．
# 仮想データ
y01 <- sample(x = 1:7, size = 1e2, replace = T)
X01 <- y01 / 1.5 + rnorm(n = 1e2, sd = 1)
plot(x = X01, y = y01)
# 上記関数による推定
irls_poisson(y_poi = y01, X_poi = X01)## $coefficient
## coef SE z
## [1,] 0.5927572 0.11681382 5.074375
## [2,] 0.2563512 0.03124537 8.204454
##
## $iteration
## [1] 10# 答え合わせ
summary(glm(formula = y01 ~ X01, family = poisson))##
## Call:
## glm(formula = y01 ~ X01, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.99939 -0.58314 0.00441 0.41091 1.63709
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.59276 0.11681 5.074 3.89e-07 ***
## X01 0.25635 0.03125 8.204 2.32e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 123.134 on 99 degrees of freedom
## Residual deviance: 55.164 on 98 degrees of freedom
## AIC: 367.08
##
## Number of Fisher Scoring iterations: 4Negative Binomial Regression
過分散のケースなどポアソン分布の仮定が妥当でない場合は負の二項分布（negative binomial distribution）へのあてはめを試してみる．コーディングに際して Silva and Rodrigues (2014, Statistics and Computing) を参考にしたとの由．
irls_nb <- function(y_nb, X_nb, maxit = 500, tol = 1e-08, offset = 0) {
# maxit （収束しない場合の）iteration回数
# tol 収束判定に用いるtolerance
y_nb <- as.vector(y_nb) # 被説明変数
X_nb <- as.matrix(X_nb) # 説明変数
reg.X <- cbind(1, X_nb) # +定数項
beta_nb <- rep(0, ncol(reg.X)) # 回帰係数の初期値
r <- 1 # 過分散パラメータ (theta) の初期値
# 「givenなrのもとで回帰係数をIRLS法で求める → givenな平均項muの理論値のもとで過分散パラメータrをNR法で求める」 を実行する繰り返し計算
for(i in 1:maxit){
for(j in 1:maxit){ # givenのrのもとで回帰係数をIRLS法で求める為の繰り返し計算
eta <- reg.X %*% beta_nb + offset # 線形予測子
mu <- exp(eta) # 平均項
w <- 1/((mu + (1/r)*mu*mu) * (1/mu)^2) # 重み行列の対角項を格納したベクトル
W <- diag(nrow(reg.X)) # n次単位行列
diag(W) <- w # n次単位行列の対角項にwで定義したベクトルの要素を代入し重み行列完成
z <- eta + (y_nb - mu)*(1/mu) - offset # 調整被説明変数(adjusted dependent variable)
old.beta <- beta_nb
beta_nb <- solve(crossprod(reg.X, W %*% reg.X)) %*% crossprod(reg.X, W %*% z) # 回帰係数を更新
if(sqrt(crossprod(beta_nb - old.beta)) < tol) break # 収束したら終了
}
U <- sum(digamma(r+y_nb) - digamma(r) + log(r) + 1 - log(r+mu) - (r+y_nb) / (r+mu)) # 対数尤度関数を過分散パラメータについて1階微分したものに前行までの繰り返し計算で得られた mu を代入
H <- sum(trigamma(r+y_nb) - trigamma(r) + (1/r) - (2/(r+mu)) + (r+y_nb) / (r+mu)^2) # 〃 r について2階微分したもの（ヘッシアン）に 〃
old.r <- r
r <- r - U/H # 過分散パラメータを更新
if (abs(old.r - r) < tol) break # 収束したら終了
}
beta_nb <- as.vector(beta_nb) # 回帰係数
var.beta <- solve(crossprod(reg.X, W %*% reg.X)) # 回帰係数の分散
stderr <- sqrt(diag(var.beta)) # 標準誤差
z.value <- beta_nb / stderr # z値
coef <- cbind(coef = beta_nb, SE = stderr, z = z.value)
s.e._of_dispersion <- sqrt(-1/H) # 過分散パラメータの標準誤差
list(coefficient = coef, iteration = i, dispersion = as.numeric(r), SE = s.e._of_dispersion) # 結果の出力
}デフォルトで用意されているデータセットMASS::quine（学校の欠席日数データ）で計算してみる．
library(MASS)
X_quine <- 1 * cbind(EthN = quine$Eth == "N",
AgeF1 = quine$Age == "F1", AgeF2 = quine$Age == "F2", AgeF3 = quine$Age == "F3")
# 上記関数による推定
irls_nb(y_nb = quine$Days, X_nb = X_quine)## $coefficient
## coef SE z
## 3.0382232 0.1956522 15.5286970
## EthN -0.5610875 0.1546672 -3.6277077
## AgeF1 -0.3855325 0.2274126 -1.6952998
## AgeF2 0.1846080 0.2312754 0.7982175
## AgeF3 0.2550005 0.2406598 1.0595891
##
## $iteration
## [1] 5
##
## $dispersion
## [1] 1.249143
##
## $SE
## [1] 0.1568657# 答え合わせ
summary(glm.nb(Days ~ Eth + Age, data = quine))##
## Call:
## glm.nb(formula = Days ~ Eth + Age, data = quine, init.theta = 1.249142793,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7902 -0.9310 -0.3795 0.3845 2.6038
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0382 0.1957 15.529 < 2e-16 ***
## EthN -0.5611 0.1547 -3.628 0.000286 ***
## AgeF1 -0.3855 0.2274 -1.695 0.090019 .
## AgeF2 0.1846 0.2313 0.798 0.424744
## AgeF3 0.2550 0.2407 1.060 0.289332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.2491) family taken to be 1)
##
## Null deviance: 192.04 on 145 degrees of freedom
## Residual deviance: 167.84 on 141 degrees of freedom
## AIC: 1107.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.249
## Std. Err.: 0.157
##
## 2 x log-likelihood: -1095.801Spatial econometrics model
SEM
\[ \ln L = -\frac{n}{2} \log(\sigma_{ML}) + \log(\det(I - \lambda W)) - \frac{n}{2} \log(2 \pi) - \frac{n}{2} \]
\[ B = I - \lambda W, \beta = \frac{ (B X)' B y }{ (B X)' B X } \]
sem_alt <- function (y_sem, X_sem, W_sem, maxit = 10000, tol = 1e-8) {
# y_sem 被説明変数ベクトル
# X_sem 説明変数行列
# W_sem 空間重み行列
# maxit （収束しない場合の）iterationの最大回数
# tol 収束判定に利用する計算誤差のtorelance
n <- nrow(X_sem) # サンプルサイズ
resid_g <- resid(lm(y_sem ~ X_sem)) # 回帰残差（初期値）
I_sem <- diag(n) # 単位行列
reg.X <- cbind(1, X_sem) # + 定数項
LL <- function (lambda_LL){ # 対数尤度関数
B <- I_sem - lambda_LL * W_sem
sigma_ML <- crossprod(B %*% resid_g) / n
objFCN <- -(n/2) * log(sigma_ML) + log(det(B)) - ((n/2) * log(2 * pi)) - n/2
objFCN
}
for (i in 1:maxit) {
optLL <- optimize(f = LL, lower = -0.99, upper = 0.99, maximum = TRUE)
lambda_ml <- optLL$maximum # 最大尤度におけるλ
l_cons <- optLL$objective # 最大化された対数尤度 max ln(L)
old.resid_g <- resid_g
IlW_sem <- I_sem - lambda_ml * W_sem # (I - λW)
beta_g <- solve(crossprod(IlW_sem %*% reg.X)) %*% t(IlW_sem %*% reg.X) %*% IlW_sem %*% y_sem
resid_g <- y_sem - reg.X %*% beta_g # 回帰残差の更新
if (sqrt(crossprod(old.resid_g - resid_g)) < tol) break # 収束したら終了
}
beta_ml <- as.numeric(beta_g)
names(beta_ml) <- c("(Intercept)", colnames(X_sem))
RSS_ml <- as.vector(crossprod(y_sem - reg.X %*% beta_ml)) # 残差二乗和
sq_sigma_ml <- RSS_ml / n # 誤差分散の最尤推定量
var_beta <- sq_sigma_ml * solve(crossprod(reg.X))
SE_ml <- sqrt(diag(var_beta)) # 回帰係数の標準誤差
t_ml <- beta_ml / SE_ml # t値
coef_ml <- cbind(coef = beta_ml, SE = SE_ml, t = t_ml)
AIC_ml <- -2 * l_cons + 2 * (ncol(reg.X)+2)
# 推定結果の出力
list(coefficient = coef_ml, lambda = lambda_ml, LL = l_cons, AIC = AIC_ml, iteration = i)
}計算してみる．
library(spdep)## Loading required package: sp## Loading required package: Matrixdata(columbus, package = "spdep")
# 重み行列の作成
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) # 行基準化
# 説明変数行列
x_col <- cbind(INC = columbus$INC, HOVAL = columbus$HOVAL)
# 推定
sem_alt_result <- sem_alt(y_sem = columbus$CRIME, X_sem = x_col, W_sem = w_std_col)
sem_alt_result## $coefficient
## coef SE t
## (Intercept) 54.6638829 5.0548580 10.814128
## INC -0.8706943 0.3566653 -2.441208
## HOVAL -0.2560451 0.1101586 -2.324331
##
## $lambda
## [1] 0.6563154
##
## $LL
## [,1]
## [1,] -181.2195
##
## $AIC
## [,1]
## [1,] 372.4391
##
## $iteration
## [1] 18SLM
\[ \ln L = -\frac{n}{2} \log(\sigma_{ML}) + \log(\det(I - \rho W)) - \frac{n}{2} \log(2 \pi) - \frac{n}{2} \]
slm_alt <- function(y_slm, X_slm, W_slm) {
n <- nrow(X_slm) # サンプルサイズ
resid_0 <- resid(lm(y_slm ~ X_slm)) # 回帰残差（初期値）
resid_L <- resid(lm(W_slm %*% y_slm ~ X_slm))
I_slm <- diag(n) # 単位行列
LL <- function (rho_LL){ # 対数尤度関数
sigma_ML <- (crossprod(resid_0 - rho_LL * resid_L)) / n
objFCN <- -(n/2) * log(sigma_ML) + log(det(I_slm - rho_LL * W_slm)) - ((n/2) * log(2 * pi)) - n/2
objFCN
}
optLL <- optimize(f = LL, lower = -0.99, upper = 0.99, maximum = TRUE)
rho_ml <- optLL$maximum # 最大尤度におけるρ
l_cons <- optLL$objective # 最大化された対数尤度 max ln(L)
lm_lag <- lm(((I_slm - rho_ml * W_slm) %*% y_slm) ~ X_slm) # ここが推定のメイン
beta_lm_lag <- coef(lm_lag)
names(beta_lm_lag) <- c("(Intercept)", colnames(X_slm))
RSS_ml <- deviance(lm_lag)
sq_sigma_ml <- RSS_ml / n # 誤差分散の最尤推定量
X_const <- cbind(1, X_slm)
var_beta <- sq_sigma_ml * solve(crossprod(X_const))
SE_ml <- sqrt(diag(var_beta)) # 回帰係数の標準誤差
t_ml <- beta_lm_lag / SE_ml # t値
coef_ml <- cbind(coef = beta_lm_lag, SE = SE_ml, t = t_ml)
AIC_ml <- -2 * l_cons + 2 * (ncol(X_const)+2)
# 推定結果の出力
list(coefficient = coef_ml, rho = rho_ml, LL = l_cons, AIC = AIC_ml)
}計算してみる．
slm_alt_result <- slm_alt(y_slm = columbus$CRIME, X_slm = x_col, W_slm = w_std_col)
slm_alt_result## $coefficient
## coef SE t
## (Intercept) 38.5515606 3.8141295 10.107565
## INC -0.8818812 0.2691208 -3.276897
## HOVAL -0.2605355 0.0831199 -3.134454
##
## $rho
## [1] 0.5154635
##
## $LL
## [,1]
## [1,] -179.881
##
## $AIC
## [,1]
## [1,] 369.7619（以上）