MCMCpackで順序ロジットのベイズ推定をするかなと思ったのですが、用意されていたのは順序プロビットのMCMCoprobitだったので、順序プロビットをします。順序ロジットの尤度関数はそんなに複雑ではないのですが、微積分をしない限りは大差ないので。
限界効果を計算する関数
分布を扱う関係でループが多く煩雑になっていますが、以下となります。
listDummies <- function(data, ndm) colnames(ndm)[!colnames(ndm) %in% colnames(data)][-1] listVariables <- function(vns, ndm, dummies){ if(is.null(vns)){ # デフォルトは全変数 vns <- colnames(ndm)[-1] } else { tmp <- c() for(vn in vns){ # データフレームを展開した行列の列名にあればそのまま使う if(vn %in% colnames(ndm)){ tmp <- c(tmp, vn) } else { # 無ければダミー変数として、展開された変数名を入れる len <- length(vn) tmp <- c(tmp, dummies[substring(dummies, 1, len) == vn]) } } vns <- tmp } # クロス項を除外しておく vns[grep("^[^:]+:[^:]+$", vns, invert = TRUE)] } getvars <- function(frml){ avars <- all.vars(frml) lst <- as.list(frml) evars <- all.vars(parse(text = lst[[length(lst)]])) dvars <- avars[!avars %in% evars] list(all=avars, explanatory=evars, dependent=dvars) } df2matrix <- function(model, data, newdata){ if(!is.null(newdata)){ vars <- getvars(model) # データフレームに余分な列があったら除外する data <- data[, vars$all] newdata <- newdata[, colnames(newdata) %in% vars$all] # データフレームに足りない列があったらエラー lacks <- vars$explanatory[!vars$explanatory %in% colnames(newdata)] if(0 < length(lacks)){ stop(paste("newdataの変数が不足しています: ", paste(lacks, collapse = ","), sep = ":")) } # もしくは足す lacks <- vars$dependent[!vars$dependent %in% colnames(newdata)] for(i in lacks) newdata[, i] <- data[1, i] # 推定にも予測にも使わないので、適当な値を入れておく mdata <- rbind(data, newdata) } else { mdata <- data } # モデル式からデータセットをつくることで、character型をダミーに展開 # 結果オブジェクトからダミー変数の取りうる値が分からないので、全データをなめる mcd <- model.matrix(model, data = mdata) if(is.null(newdata)) mcd <- colMeans(mcd) else mcd <- mcd[1:nrow(newdata), ] if(!is.matrix(mcd)) mcd <- matrix(mcd, 1, dimnames = list(NULL, names(mcd))) mcd } predict.MCMCoprobit <- function(model, data, posterior, newdata, IsBetaMeans = TRUE){ if("mcmc" != class(posterior)) stop("MCMCoprobitの推定結果でないと処理できません．") # newdataが行列であれば、df2matrixで処理済みと見なす if(is.matrix(newdata)) mcd <- newdata else { if(is.null(data)) stop("dataがnullです．") mcd <- df2matrix(model, data, newdata) } # 推定結果の係数 if(IsBetaMeans) coef <- colMeans(posterior) else coef <- posterior if(!is.matrix(coef)) coef <- matrix(coef, ncol = 1, dimnames = list(names(coef), NULL)) else coef <- t(coef) coef <- rbind(coef, rep(0, rep(ncol(coef)))) rownames(coef)[nrow(coef)] <- "gamma1" betas <- colnames(mcd) gammas <- rownames(coef)[!rownames(coef) %in% betas] # (線形結合, 係数) Xb <- mcd %*% as.matrix(coef[betas, ]) # (観測値（i.e. dataのレコード番号）, 係数（i.e. MCMCで推定したもの）, カテゴリー)ごとの確率 lapply(1:nrow(Xb), function(i){ p <- sapply(sort(gammas), function(j) pnorm(coef[j, ] - Xb[i, ]), USE.NAMES = FALSE) if(!is.matrix(p)) p <- matrix(p, 1) p <- cbind(rep(0, nrow(p)), p, rep(1, nrow(p))) q <- sapply(2:ncol(p), function(k) p[, k] - p[, k - 1]) if(!is.matrix(q)) return(matrix(q, 1)) q }) } calcMarginalEffect <- function(model, posterior, data, newdata = NULL, vns = NULL, h = 1e-5, a = 0.025, level = 0.95, rev.dum = TRUE, levels = NULL){ # 限界効果を計測する点 ndm <- df2matrix(model, data, newdata) # ダミー変数を特定する dummies <- listDummies(data, ndm) # 説明変数の行列のある列の値を変える # 単項だけではなく、クロス項（e.g. x*z）の部分も変化させる terms <- strsplit(colnames(ndm), ":") chgvalue <- function(m, name, value){ m2 <- m for(i in 1:ncol(m)){ cns <- terms[[i]] r <- rep(1, nrow(m)) for(j in 1:length(cns)){ r <- r*((cns[j] != name) * m[, cns[j]] + (cns[j] == name) * value) } m2[, i] <- r } m2 } vns <- listVariables(vns, ndm, dummies) # 変数名→{newdata}の行番号→(係数, 応答)行列 me <- lapply(vns, function(vn){ if(rev.dum && vn %in% dummies){ ndm2 <- rbind(chgvalue(ndm, vn, 1), chgvalue(ndm, vn, 0)) h <- 0.5 } else { ndm2 <- rbind(chgvalue(ndm, vn, ndm[, vn] + h), chgvalue(ndm, vn, ndm[, vn] - h)) } prdct <- predict.MCMCoprobit(model, NULL, posterior, ndm2, IsBetaMeans = FALSE) r <- lapply(1:nrow(ndm), function(i){ r <- (prdct[[i]] - prdct[[i + nrow(ndm)]])/(2*h) colnames(r) <- levels r }) r }) names(me) <- vns a <- (1 - level)/2 qntl <- lapply(vns, function(vn){ r <- lapply(me[[vn]], function(m){ r <- t(sapply(1:ncol(m), function(j){ quantile(m[, j], c(0, a, 0.5, (1-a), 1)) })) rownames(r) <- levels r }) r }) names(qntl) <- vns means <- lapply(vns, function(vn){ r <- t(sapply(me[[vn]], function(m){ colMeans(m) })) }) names(means) <- vns list(mean = means, quantile = qntl, each = me) }
データセット
順序ロジットで使ったものと同じデータ生成プロセスを使います。
df01 <- with(new.env(), { n <- 100 x1 <- runif(n) x2 <- x1/2 + runif(n) d <- letters[round(runif(n, min = 0.5, max = 3.5))] y <- 1 + 2*x1 - 3*x2 + (d == "b")*1 + rnorm(n) g <- as.ordered(0 + (y >= 0)*1 + (y >= 1)*1 + (y >= 2)*1) data.frame(g, x1, x2, d) })
順序プロビットのベイズ推定
水準あたりは無くてもよいのですが、結果が分かりやすいので。
library(MCMCpack) # 名前をつけるのに使うので、水準を保存しておく levels <- sprintf("response: %s", levels(df01$g)) # MCMCoprobitがordered型だと警告を出すので数値型に変える o2n <- function(a) (1:length(levels(a)))[a] df01$g <- o2n(df01$g) # 推定モデルは変数にとっておく model <- g ~ x1 + x2 + d # 推定する閾値の数を求めておく n_a0 <- length(levels) - 1 # 薄く広い事前分布を設定して推定 # b0, B0, a0, A0は無くても推定できるが、ベイズファクターが正しく使えなくなる r_mcmc <- MCMCoprobit(model, data = df01, b0 = c(1, 1, -1, 1, -1), B0 = rep(1/10, 5), a0 = 1:n_a0, A0 = rep(1/10, n_a0)) # 推定結果の要約 s <- summary(r_mcmc)
限界効果の計算
MCMCoprobitの推定結果（事後分布）ができていれば、上の関数に投げれば済みます。
# 個体の値が平均値だとした限界効果の計算 # 推定に使ったmodelとデータも必要なので注意 me <- calcMarginalEffect(model, r_mcmc, df01, levels = levels)
結果が分布で出てくるのでややこしいのですが、平均にしてしまえば話がしやすくなります。説明変数ごとにカテゴリーの確率ごとの限界効果が出ます。
me$mean
分布があるので95%信用区間（というか変則分位数）の計算も同時に行っています。
me$quantile$x1
平均限界効果の計算
データセットの行数×サンプリングの行数の膨大な事後分布ができて計算機リソースを消耗するのでお勧めしませんが、原理的には簡単にできます。 # 行番号17から18の変数x1とx2の平均限界効果を求める me <- calcMarginalEffect(model, r_mcmc, df01, vns = c("x1", "x2"), newdata = df01[17:19, ], levels = levels) lapply(me$each, function(v){ r <- v[[1]] for(k in 2:length(v)){ r <- rbind(r, v[[k]]) } colMeans(r) })
ご利用上の注意
検算をしていないので、誰か代わりに('-' )＼(--；)BAKI
