Data Analysis Using Statistical Packages: Discrete Regression Analysis
Sho Kuroda / 黒田翔
2024年3月 (Last update: 2024年4月)
library(tidyverse)1 Maximum likelihood / 最尤法
1.1 Bernoulli distribution / ベルヌーイ分布
5回の試行で3回成功（成功確率 \(p\)），2回失敗（失敗確率 \(1-p\)）した場合の尤度関数と対数尤度関数．
\[ L(p) = \Pr(\mbox{failure})^2 \times \Pr(\mbox{success})^3 = (1-p)^2 p^3, \quad \log L(p) = 2 \log (1-p) + 3 \log (p) \]
一階の条件（\(L(\hat{p})' = 0\) または \(\log L(\hat{p})' = 0\)）より \(\hat{p} = 0.6\) が得られる．
\[ \frac{d L(p)}{d p} = \frac{d}{d p} (1-p)^2 p^3 = -2(1-p)p^3 + (1-p)^2 \times 3p^2 = (1-p)p^2[-2p + 3(1-p)] = 0 \]
l_func <- function (p) (1-p)^2 * p^3 # 尤度関数
ll_func <- function (p) 2 * log(1-p) +3 * log(p) # 対数尤度関数
plot(l_func, xlim = c(0, 1)); abline(v = .6, lty = 2, col = 2)
plot(ll_func, xlim = c(0, 1)); abline(v = .6, lty = 2, col = 2)
1.2 Normal distribution / 正規分布
標本 \(\{-3, -2, 0.5, 0.5, 1.5\}\) に基づいて母平均 \(\mu\) と母分散 \(\sigma^2\) を推定．
この標本が正規母集団からサンプリングされていると仮定すると，対数尤度関数は次のように書ける．
\[ f = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left[- \frac{(y_i - \mu)^2}{2 \sigma^2} \right] , \quad \log f = - \frac{1}{2} \log 2 \pi - \frac{1}{2} \log \sigma^2 - \frac{(y_i - \mu)^2}{2 \sigma^2} \]
\[ \log L = - \frac{n}{2} \log 2 \pi - \frac{n}{2} \log \sigma^2 - \frac{1}{2 \sigma^2} \sum_i (y_i - \mu)^2 \]
よって，最尤推定量は次のように解析的に求められる．
\[ \hat{\mu} = \frac{1}{n} \sum_i y_i , \quad \hat{\sigma}^2 = \frac{1}{n} \sum_i (y_i - \bar{y})^2 \]
y <- c(-3, -2, 0.5, 0.5, 1.5) # 標本
mean(y) # 標本平均（mu の最尤推定量と同じ）## [1] -0.5var(y) # 標本分散## [1] 3.625sum((y - mean(y))^2) / length(y) # 分散（sigma^2 の最尤推定量と同じ）## [1] 2.9対数尤度関数を定義する． 最初の引数がパラメタベクトルになっている（後述の最適化のため）．
ll_func_n <- function (param, y, n) { # 対数尤度関数
mu <- param[1]
sigma2 <- param[2]
-n/2*log(2*pi) - n/2*log(sigma2) - 1/(2*sigma2)*sum((y-mu)^2)
}\(\sigma^2 = 1\) に固定した場合の \(\mu\) と対数尤度の関係をプロットする．
mu_list <- seq(-2, 2, length = 100)
ll_mu_list <- sapply(X = mu_list,
FUN = function(x) ll_func_n(param = c(x, 1), y = y, n = length(y)))
plot(x = mu_list, y = ll_mu_list, xlim = c(-2, 2), type = "l")
3Dで描画する．（outer 関数を使って z
を計算するつもりだったが，なぜか outer
関数にかませる尤度関数中で sum((y-mu)^2) が評価できないので
for 文をネストして計算する．）
grid_mu <- seq(-2, 2, length = 20) # 探索する範囲を grid で表す
grid_sigma2 <- seq(1, 10, length = 20)
z <- matrix(0, nrow = length(grid_mu), ncol = length(grid_sigma2))
for (i in 1:length(grid_mu)) { # 探索 grid に対応する対数尤度を計算
for (j in 1:length(grid_sigma2)) {
z[i, j] <- ll_func_n(param = c(grid_mu[i], grid_sigma2[j]), y = y, n = length(y))
}
}
which_max_ml <- which(z == max(z), arr.ind = T) # 尤度を最大化する grid の index
grid_mu[which_max_ml[1]]; grid_sigma2[which_max_ml[2]] # ML推定量 （gridが粗いので精度は低い）## [1] -0.5263158## [1] 2.894737persp(grid_mu, grid_sigma2, z, theta = 50, phi = 20) -> pmat_persp
points(trans3d(x = grid_mu[which_max_ml[1]], y = grid_sigma2[which_max_ml[2]],
z = z[which_max_ml], pmat = pmat_persp), col = 2, pch = 19)
汎用の最適化関数 optim を使って最尤推定する．
$par が最適化されたパラメタの値（最尤推定値）．
$value は最大化された対数尤度．
$convergence は 0 ならば収束したことを意味する．
optim
はデフォルトでは与えられた関数を最小化するため，最大化するためには目的関数を
-1 倍した値を最小化するように設定する（引数
control = list(fnscale = -1)）．
optim(par = c(0, 1), fn = ll_func_n, y = y, n = length(y), control = list(fnscale = -1))## $par
## [1] -0.5000029 2.8995297
##
## $value
## [1] -9.75647
##
## $counts
## function gradient
## 65 NA
##
## $convergence
## [1] 0
##
## $message
## NULL1.3 Linear regression / 線形回帰モデル
誤差項が正規分布に従うと仮定した場合の回帰モデルの対数尤度関数．
\[ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + u_i, \quad u_i \sim N(0, \sigma^2) \]
\[ \log L (\beta) = - \frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum_i (y_i - x_i' \beta)^2 \]
ここで \(x_{i}' \beta = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}\) である．
引数 param は k+1 個のパラメタベクトル． k
は定数項を含めた説明変数の数． 1 つ目から k 個目までが回帰係数で，k+1
個目は \(\sigma^2\)．
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])
}
n <- nrow(x_gauss) # サンプルサイズ
k <- ncol(x_gauss) # 説明変数の数（定数項を含む）
optim(par = rep(1, k+1), fn = LL_gauss, control = list(fnscale = -1)) # 尤度最大化
}
y_swiss <- swiss$Fertility # y
x_swiss <- cbind(1, swiss$Agriculture, swiss$Examination) # X
ml_gauss(y_gauss = y_swiss, x_gauss = x_swiss)## $par
## [1] 94.66923521 -0.09443892 -1.19762789 86.69205312
##
## $value
## [1] -171.5454
##
## $counts
## function gradient
## 357 NA
##
## $convergence
## [1] 0
##
## $message
## NULL答え合わせ．
lm(Fertility ~ Agriculture + Examination, swiss)##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination, data = swiss)
##
## Coefficients:
## (Intercept) Agriculture Examination
## 94.610 -0.094 -1.1952 Discrete regression model / 離散回帰モデル
Titanic データの読み込み．
titanic <- read.csv("https://raw.githubusercontent.com/kurodaecon/bs/main/data/titanic3_csv.csv")2.1 Linear probability model / 線形確率モデル
生存したかどうかを表す Survived
変数を，性別，年齢，客室等級変数に回帰．
\[ \mbox{Survived} = \beta_0 + \beta_1 \mbox{Male} + \beta_2 \mbox{Age} + \beta_3 \mbox{2nd-class} + \beta_4 \mbox{3rd-class} \]
summary(lm(survived ~ sex + age + factor(pclass), data = titanic))##
## Call:
## lm(formula = survived ~ sex + age + factor(pclass), data = titanic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09442 -0.24408 -0.08375 0.23289 0.99397
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1049549 0.0438211 25.215 < 2e-16 ***
## sexmale -0.4914132 0.0255520 -19.232 < 2e-16 ***
## age -0.0052695 0.0009316 -5.656 2.00e-08 ***
## factor(pclass)2 -0.2113738 0.0348568 -6.064 1.85e-09 ***
## factor(pclass)3 -0.3703874 0.0325039 -11.395 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3913 on 1041 degrees of freedom
## ( 263 個の観測値が欠損のため削除されました )
## Multiple R-squared: 0.3691, Adjusted R-squared: 0.3667
## F-statistic: 152.3 on 4 and 1041 DF, p-value: < 2.2e-16解釈
- 男性の生存率は49%ポイント低い
- 年齢が1歳高い乗客の生存率は0.5%ポイント低い
- 2等客室の乗客の生存率は1等と比べて21%ポイント低い
2.2 Binary probit model / 二値プロビット・モデル
\[ p = \Pr (y = 1) = \Phi (\beta_0 + \beta_1 x) \]
glm 関数で推定．
probit_titanic <- glm(formula = survived ~ sex + age + factor(pclass), data = titanic,
family = binomial(link = "probit"))
summary(probit_titanic)##
## Call:
## glm(formula = survived ~ sex + age + factor(pclass), family = binomial(link = "probit"),
## data = titanic)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.055101 0.181290 11.336 < 2e-16 ***
## sexmale -1.485640 0.093813 -15.836 < 2e-16 ***
## age -0.019426 0.003606 -5.387 7.16e-08 ***
## factor(pclass)2 -0.760170 0.129921 -5.851 4.89e-09 ***
## factor(pclass)3 -1.303160 0.126633 -10.291 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1414.62 on 1045 degrees of freedom
## Residual deviance: 984.55 on 1041 degrees of freedom
## ( 263 個の観測値が欠損のため削除されました )
## AIC: 994.55
##
## Number of Fisher Scoring iterations: 52.2.1 Estimation by numerical calculation / 数値計算による推定
対数尤度関数を汎用の最適化関数 optim
で最大化することで推定する．
選択確率は以下で与えられる．
\[ \Pr (Y=1) = \Phi (\beta_0 + \beta_1 X) = \frac{1}{2 \pi} \int_{-\infty}^{\beta_0 + \beta_1 X} \exp \left( - \frac{t^2}{2} \right) dt \]
よって，尤度関数は以下のようになる．
\[ L (\beta_0, \beta_1) = \prod_i \left[ \{ \Pr (Y_i = 1) \}^{Y_i} \times \{ \Pr (Y_i = 0) \}^{1 - Y_i} \right] = \prod_i \left[ \{ \Phi (\beta_0 + \beta_1 x) \}^{Y_i} \times \{ 1 - \Phi (\beta_0 + \beta_1 x) \}^{1 - Y_i} \right] \]
対数を取ると対数尤度関数が得られる．
\[ \log L (\beta_0, \beta_1) = \sum_i \left[ Y_i \times \log \{ \Phi (\beta_0 + \beta_1 x) \} + (1 - Y_i) \times \log \{ 1 - \Phi (\beta_0 + \beta_1 x) \} \right] \]
これを，\((\beta_0, \beta_1)\) を引数とする関数として定義すればよい．
標準正規分布の分布関数は 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))
}上と同じ説明変数を使って推定する．
まずは，計算のために character/factor 型に対応する数値を計算しておく．
titanic_numeric <- titanic %>%
filter(!is.na(age)) %>%
mutate(male = 1*(sex == "male"),
pclass2 = 1*(pclass == 2),
pclass3 = 1*(pclass == 3))titanic_x <- titanic_numeric %>%
dplyr::select(male, age, pclass2, pclass3) %>%
as.matrix
ml_probit(y_probit = titanic_numeric$survived, x_probit = cbind(1, titanic_x))## $par
## [1] 2.05418759 -1.48531017 -0.01941115 -0.76010160 -1.30246558
##
## $value
## [1] -492.2765
##
## $counts
## function gradient
## 468 NA
##
## $convergence
## [1] 0
##
## $message
## NULL2.2.2 Marginal effect / 限界効果
説明変数の平均値 \(\bar{x}\) を用いる方法．
\[ ME_{x} = \phi (\hat{\beta_0} + \hat{\beta_1} \bar{x}) \hat{\beta_1} \]
xb_mean <- probit_titanic$coef %*% c(1, mean(titanic_numeric$male), mean(titanic$age, na.rm = TRUE),
mean(titanic$pclass == 2), mean(titanic$pclass == 3))
as.numeric(dnorm(x = xb_mean)) * probit_titanic$coef## (Intercept) sexmale age factor(pclass)2 factor(pclass)3
## 0.77727888 -0.56189779 -0.00734726 -0.28751081 -0.49288013Individual ごとの限界効果の平均値を計算する方法．
\[ ME_{x} = \frac{1}{n} \sum_i ME_{i, x} = \frac{1}{n} \sum_i \phi (\hat{\beta_0} + \hat{\beta}_1 x_i) \hat{\beta}_1 \]
xb <- as.matrix(cbind(1, titanic_numeric[, c("male", "age", "pclass2", "pclass3")])) %*% probit_titanic$coef
c(male = mean(dnorm(x = xb) * probit_titanic$coef[2], na.rm = TRUE),
age = mean(dnorm(x = xb) * probit_titanic$coef[3], na.rm = TRUE),
class2 = mean(dnorm(x = xb) * probit_titanic$coef[4], na.rm = TRUE),
class3 = mean(dnorm(x = xb) * probit_titanic$coef[5], na.rm = TRUE))## male age class2 class3
## -0.392874167 -0.005137142 -0.201025116 -0.344617610mfx パッケージを使う場合．
# install.packages("mfx")
library(mfx)
probitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)## Call:
## probitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## sexmale -0.5408979 0.0287920 -18.7864 < 2.2e-16 ***
## age -0.0074648 0.0013850 -5.3898 7.052e-08 ***
## factor(pclass)2 -0.2672290 0.0404958 -6.5989 4.141e-11 ***
## factor(pclass)3 -0.4666238 0.0395420 -11.8007 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "sexmale" "factor(pclass)2" "factor(pclass)3"2.3 Binary logit model / 二値ロジット・モデル
\[ p = \Pr (y = 1) = \Lambda (\beta_0 + \beta_1 x) \]
glm 関数で推定．
logit_titanic <- glm(formula = survived ~ sex + age + factor(pclass), data = titanic,
family = binomial(link = "logit"))
summary(logit_titanic)##
## Call:
## glm(formula = survived ~ sex + age + factor(pclass), family = binomial(link = "logit"),
## data = titanic)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.522075 0.326702 10.781 < 2e-16 ***
## sexmale -2.497845 0.166037 -15.044 < 2e-16 ***
## age -0.034393 0.006331 -5.433 5.56e-08 ***
## factor(pclass)2 -1.280571 0.225538 -5.678 1.36e-08 ***
## factor(pclass)3 -2.289661 0.225802 -10.140 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1414.62 on 1045 degrees of freedom
## Residual deviance: 982.45 on 1041 degrees of freedom
## ( 263 個の観測値が欠損のため削除されました )
## AIC: 992.45
##
## Number of Fisher Scoring iterations: 4余裕がある履修者向けの宿題：プロビットと同様に対数尤度関数を書いて
optim で最適化計算する．
2.3.1 Marginal effect / 限界効果
logitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)## Call:
## logitmfx(formula = survived ~ sex + age + factor(pclass), data = titanic)
##
## Marginal Effects:
## dF/dx Std. Err. z P>|z|
## sexmale -0.5514348 0.0293172 -18.8092 < 2.2e-16 ***
## age -0.0080960 0.0014881 -5.4405 5.312e-08 ***
## factor(pclass)2 -0.2673472 0.0404682 -6.6064 3.939e-11 ***
## factor(pclass)3 -0.4901785 0.0400797 -12.2301 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "sexmale" "factor(pclass)2" "factor(pclass)3"余裕がある履修者向けの宿題：プロビットと同様に（mfx::logitmfx
関数を使わずに）計算する．
2.4 Advanced: Conditional logit
演習では扱わないが，興味がある履修者は 自動車保有の離散選択問題：Conditional logit を試してみよう．
3つ以上の選択肢から1つを選ぶ問題も上記を一般化して定式化・推定できることが分かるだろう．
.