アドベントカレンダー
株式会社GENDAでデータサイエンティストをしているtoma2です。
この記事は、GENDAアドベントカレンダー2024の3日目の記事になります。
GENDAアドベントカレンダーでは、プロダクトや組織開発に関わるメンバーを中心に多様なテーマの記事を投稿しています。ぜひ、購読登録をしていただき12月25日までお楽しみください。
はじめに
最近、一般線形モデルの予測値を合計することがありました。
その際、予測値の合計の予測区間を求めたかったのですが、RやPythonにそのような関数はなく、日本語での解説記事もWeb上では見当たりませんでした。
Stack Overflowを探したところ、以下の質問で解説されており、数式での説明と実用的な実装が紹介されていました。回帰分析の基礎がわかっていれば、簡単に計算できるものでした。
この記事では、Stack Overflowの回答を参考にダミーデータを用いて、予測値の合計の予測区間の算出方法を紹介します。なお、回答で実用的な実装で計算が追いにくくなっている箇所は、回帰分析の復習も兼ねて愚直な実装にしました。最後にベイズ流回帰での結果と比較します。
ダミーのデータとシチュエーション
身長と体重のダミーデータを利用して計算していきます。
90人分の身長と体重のデータ（学習データ）から回帰モデルを作成し、身長のみの10人のデータ（テストデータ）の合計体重の予測値と予測区間を求めるようなシチュエーションを想定してみます。
library(MASS)
# パラメータの設定
# 身長と体重のパラメータはe-stat(https://www.e-stat.go.jp/dbview?sid=0003224177)から引用
set.seed(123)
n <- 100
height_mu <- 167.7
height_s <- 6.9
weight_mu <- 67.4
weight_s <- 12.0
cor_hw <- 0.5 #相関は適当に0.5
# 多変量正規分布からの乱数生成
mu <- c(height_mu, weight_mu)
sigma <- matrix(c(height_s**2, cor_hw*height_s*weight_s,
cor_hw*height_s*weight_s, weight_s**2), nrow = 2)
samples <- mvrnorm(n = n, mu = mu, Sigma = sigma)
data <- data.frame(height = samples[,1], weight = samples[,2])
# 学習データと検証データに分割
train <- data[1:90, ]
test <- data[(90+1):n, 'height', drop = FALSE]
学習データのプロット
回帰モデルの作成と予測値の合計
計算は以下のように簡単にできます。
# 身長から体重への回帰
lm_model <- lm(weight~height,data=train)
summary(lm_model)
# Call:
# lm(formula = weight ~ height, data = train)
#
# Residuals:
# Min 1Q Median 3Q Max
# -20.7391 -7.2235 0.4567 7.0219 19.6649
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -62.5271 25.5109 -2.451 0.0162 *
# height 0.7741 0.1514 5.114 1.83e-06 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 9.405 on 88 degrees of freedom
# Multiple R-squared: 0.2291, Adjusted R-squared: 0.2203
# F-statistic: 26.15 on 1 and 88 DF, p-value: 1.825e-06
# 予測値の合計の計算
pred <- predict(lm_model, newdata = test)
pred_sum <- sum(pred)
print(pred_sum)
# [1] 698.761
予測値の合計の予測区間
予測値の合計の分散を計算することで、予測区間を求められます。
Stack Overflowに
Xp %*% vcov(lmObject) % t(Xp), which is slow
と記載してある愚直な方法で計算してみます。
# 残差分散
residual_var <- summary(lm_model)$sigma**2
# 予測値の分散共分散行列を求める
# vcov(lm_model)は回帰係数の分散共分散行列
xp <- cbind(1, test$height)
y_hat_var <- xp %*% vcov(lm_model) %*% t(xp)
# 予測値の合計の分散は、予測値の分散共分散行列の合計に予測値分の残差分散を足したもの
sum_var <- sum(y_hat_var) + residual_var * nrow(y_hat_var)
alpha <- 0.95
qt_value <- qt((1 - alpha) / 2, lm_model$df.residual, lower.tail = FALSE)
# 予測値の合計の予測区間
lm_result <- pred_sum + c(-1,0,1) * qt_value * sqrt(sum_var)
print(lm_result)
# [1] 635.9544 698.7610 761.5676
予測値の合計の予測区間を求めることができました。
予測値の合計の予測区間（ベイズ流回帰）
ベイズ流回帰では、予測値の合計の予測区間を簡単に求められます。検算も兼ねて計算してみます。
# ベイズ流回帰での予測値の合計の予測区間
library(rstanarm)
stan_model <- stan_glm(weight ~ height, data = train, family = gaussian(), seed = 123)
summary(stan_model)
# Model Info:
# function: stan_glm
# family: gaussian [identity]
# formula: weight ~ height
# algorithm: sampling
# sample: 4000 (posterior sample size)
# priors: see help('prior_summary')
# observations: 90
# predictors: 2
#
# Estimates:
# mean sd 10% 50% 90%
# (Intercept) -62.3 25.5 -94.9 -62.3 -29.9
# height 0.8 0.2 0.6 0.8 1.0
# sigma 9.5 0.7 8.6 9.5 10.4
#
# Fit Diagnostics:
# mean sd 10% 50% 90%
# mean_PPD 67.8 1.4 66.0 67.8 69.6
#
# The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
#
# MCMC diagnostics
# mcse Rhat n_eff
# (Intercept) 0.4 1.0 3732
# height 0.0 1.0 3710
# sigma 0.0 1.0 3430
# mean_PPD 0.0 1.0 3325
# log-posterior 0.0 1.0 1732
#
# For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
stan_pred <- posterior_predict(stan_model, newdata = test)
stan_result <- quantile(rowSums(stan_pred), probs = c(0.025, 0.5, 0.975))
print(stan_result)
# 2.5% 50% 97.5%
# 634.2317 699.1502 762.3585
もちろんlmでの結果と完全に一致はしませんが、おおよそ同じ値となっていることが確認できました。
宣伝
GENDAデータチームでは、プロダクトのデータ分析や機械学習プロジェクトを推進できるデータサイエンティストを現在募集しています。カジュアル面談も可能ですので、お気軽にご連絡ください。
Discussion