打点と得点の合計値を他の指標で回帰するベイズ推定

野球の出塁率と長打率を足し合わせた値OPSの評価をするために、打点と得点の合計値を他の指標で回帰するコードです。非負制約をつけるために、ベイズ推定をしています。

なお、結果は以下で、OLSで推定してもほとんど同様の値になりました。線形回帰モデルなんで当たり前なんですが。

	単打	二塁打	三塁打	本塁打	四死球
限界効果	0.49	0.86	1.36	2.79	0.32
単打比	1	1.74	2.75	5.66	0.66

打点と得点の合計値への限界効果

library(MCMCpack)

df01 <- read.table("https://wh.anlyznews.com/R/dataset/npb_b_ranking.txt", sep ="\t", header = TRUE)
df01$単打 <- with(df01, 安打 - 二塁打 - 三塁打 - 本塁打)

# モデル式から説明変数の行列と、被説明変数をつくる
frml <- 打点 + 得点 ~ I(単打 + 二塁打 + 三塁打 + 本塁打) +  I(二塁打 + 三塁打 + 本塁打) +  I(三塁打 + 本塁打) + 本塁打 + I(四球 + 死球) + 盗塁 + 盗塁刺 + I(犠打 + 犠飛)
df02 <- model.frame(frml, df01)
X <- model.matrix(frml, df02)
y <- model.response(df02)

# 切片項と盗塁と盗塁刺以外は、非負制約をつける（事前分布とL-BFGS-Bの下限を変える）
no.limit <- which(colnames(X) %in% c("(Intercept)", "盗塁", "盗塁刺"))

# MCMCを実行するときの初期値を決める最尤法（L-BFGS-B）で使う下限を設定
no.limit.flag <- rep(TRUE, ncol(X))
no.limit.flag[no.limit] <- FALSE
lower <- ifelse(no.limit.flag, 1e-12, Inf)
upper <- c(rep(Inf, ncol(X)))

# 初期値を設定
theta.init <- ifelse(no.limit.flag, 0.5, 0.0)

# 尤度関数
llf_nd <- function(theta){
# concentrated loglikelihood functionにして分散の推定は省略
  e <- y - X %*% theta
  sd <- sqrt(sum((e - mean(e))^2)/length(e))
  sum(dnorm(e, mean=0, sd=sd, log=TRUE))
}

# 事前分布
prior.u <- function(theta){
  # 非負制約をつけるために、事前分布は一様分布
  a <- 10
  sum(dunif(theta[no.limit], min=-a, max=a, log=TRUE) + 
	dunif(theta[-no.limit], min=0, max=a, log=TRUE))
}

# MCMCを実行
posterior <- MCMCmetrop1R(function(theta){
  v <- llf_nd(theta) + prior.u(theta)
  # -Infを返すとエラーが出るので絶対値が大きな負の数を戻す（大きすぎても特異になってエラー）
  if(!is.finite(v)){
    return(-1e+8)
  }
  return(v)
}, theta.init = rep(0.5, ncol(X)), 
	mcmc = 100000, 
	optim.method="L-BFGS-B", optim.lower=lower, optim.upper=upper)

# 推定結果のまとめの列名を書き換える
(s <- summary(posterior))
rownames(s$statistics) <- colnames(X)

# サンプリング結果から限界効果の区間推定を計算する
sp <- 2
cn <- c("単打", "二塁打", "三塁打", "本塁打", "四死球")
me <- posterior[, sp]
me <- cbind(me, apply(posterior[, sp:(sp + 1)], 1, sum))
me <- cbind(me, apply(posterior[, sp:(sp + 2)], 1, sum))
me <- cbind(me, apply(posterior[, sp:(sp + 3)], 1, sum))
me <- cbind(me, posterior[, sp + 4])
me_mean <- apply(me, 2, mean)
me_mean <- matrix(c(me_mean, me_mean/me_mean[1]), nrow = 2, byrow = TRUE, dimnames = list(c("限界効果", "単打比"), cn))
alpha <- 0.05
me_ci <- apply(me, 2, \(x) quantile(x, probs = c(0, alpha/2, 1-alpha/2, 1)))
colnames(me_ci) <- cn

write.table(me_ci, "me_ci.txt", sep = "\t")
write.table(me_mean, "me_mean.txt", sep = "\t")