前のエントリーを書いているときに、Rのordinalパッケージに限界効果を計算する関数が用意されていないことに気づきました。MASSパッケージとererパッケージを使えばよいわけですが*1、ererパッケージの更新が微妙な感じがしたので、ordinalパッケージの結果オブジェクトから限界効果を計算する関数を用意します。
順序ロジスティック回帰
順序ロジットは、①大好き、②好き、③どちらでもない、④嫌い…と言うような、順序は分かるが、数量ではないデータを扱う回帰分析です。それぞれ選択肢と言うかカテゴリーがとる確率は非負で、また合計は1になり、さらに選択肢と選択肢の差は一定とは限らないため、普通の線形回帰は用いることができないため、アンケートデータの分析の定番手法となっています。
典型的な順序データの分析で簡素なものは、個人が、
種類の選択肢の中から
を選ぶ確率
を、
と考えます。
では
は値枠が0から1までの単調増加関数で、要するに何か累積分布関数です。計算機の性能が低い時代の名残で、手計算で用意に微分ができるロジスティック分布の確率密度関数
)がよく用いられています。ただし、
のとき0、
のときは1とします。
係数ベクトルと
の推定は最尤法を用います。尤度関数は、サンプルサイズを
、
を
のとき1、
のとき0とする単関数として、
です。間違いがないようにordinal::clmやMASS::polrなどによる推定を推奨しますが、Rのnlm関数やoptim関数あたりに対数尤度関数を与えて推定することは、そんなに難しくは無いです。
問題は推定された係数ベクトルの解釈です。一次以外の項の無い線形回帰であれば推定された係数はそのまま限界効果を表します。説明変数を増減させたときの影響はすぐ分かります。しかし、非線形関数を推定するので、限界効果は個人の状態（説明変数の具体的な値）に依存し、自明ではないです。そこで、限界効果（
は説明変数の1つ）や、観測値すべての限界効果を計算した平均限界効果を計算して解釈しようと言うことになります。ここ20年ぐらいは一般的な話なのですが、意外に限界効果の計算は標準サポートされていません。
限界効果を計算する関数
単純に数値微分して計算します。平均値以外の限界効果を計算したり、限界効果のプロットをする場合などに備えて、一部の変数のものだけを計算できるようにしたり、複数の点を同時に計算できるようにしました。またMASSパッケージのpolrの結果も処理できるようにしました。ただし、linkはlogisticとprobitにしか対応していません。
# # r_obj: MASSパッケージのpolrもしくはordinalパッケージのclmの結果オブジェクト # data: 結果オブジェクトの計算に使ったデータフレーム（デフォルトはr_objより生成） # newdata: 限界効果を計算する点のデータフレーム（デフォルトは平均値） # vns: 限界効果を計算する変数（デフォルトは説明変数すべて） # IsTreatingValuesOfDummiesAsMeans: dataを省略時、ダミーは平均を値にして計算（デフォルトはTRUE)，FALSEでは最頻値で計算 # rev.dum: ダミーの限界効果の計算の補正（デフォルトはTRUE) # interval: 信頼区間も計算するかどうか（デフォルトはFALSE) # level: 信頼水準（デフォルトは0.95） # h: テイラー展開で使う微細な値 # return: 複数の点が入力されれば限界効果をまとめて計算するので、結果はリストで戻す # calcMarginalEffect <- function(r_obj, data = data.frame(r_obj$model), newdata = NULL, vns = NULL, IsTreatingValuesOfDummiesAsMeans = TRUE, rev.dum = TRUE, interval = FALSE, level = 0.95, h = 1e-5){ if(is.null(newdata)){ newdata <- data.frame(lapply(data, function(x){ if(is.numeric(x)) return(mean(x, na.rm = FALSE)) # model.matrixがNAがあるデーターフレームの行を落とすので、使わない値を入れる if(IsTreatingValuesOfDummiesAsMeans) return(x[1]) # characterは最頻値がデフォルト freq <- table(x) names(freq)[which(freq == max(freq, na.rm = FALSE))] })) } else { # データフレームの指定があれば、IsTreatingValuesOfDummiesAsMeansは無視する IsTreatingValuesOfDummiesAsMeans <- FALSE # 被説明変数の名前を取得 frml <- terms(r_obj$model) avars <- all.vars(frml) lst <- as.list(frml) evars <- all.vars(lst[[length(lst)]]) dvars <- avars[!avars %in% evars] # 被説明変数が省略されているときは、列をつくる # model.matrixがNAがあるデーターフレームの行を落とすので、使わない値を入れる for(dv in dvars){ freq <- table(data[, dv]) newdata[, dv] <- names(freq)[which(freq == max(freq, na.rm = FALSE))] } # 説明変数が足りない場合はエラーを出して止まる if(any(!avars %in% colnames(newdata))){ stop(paste("the new data.frame doesn't have the column", paste(avars[!avars %in% colnames(newdata)], collapse = ","), sep = ":")) } } # モデル式からデータセットをつくることで、character型をダミーに展開 # 結果オブジェクトからダミー変数の取りうる値が分からないので、全データをなめる mcd <- model.matrix(terms(r_obj$model), data = rbind(newdata, data)) if(nrow(mcd) != nrow(newdata) + nrow(data)) stop("data frames are mismatching!") # データフレームになくて、データセットにある列は、展開されたダミー変数 dummies <- colnames(mcd)[!colnames(mcd) %in% colnames(data)][-1] # I(...)は除く excludeME <- function(v){ n <- grep("^I\\([^)]+\\)$", v) if(0>=length(n)) return(v) v[-n] } dummies <- excludeME(dummies) # ダミー変数に平均値というか期待値を入れる if(IsTreatingValuesOfDummiesAsMeans){ k <- 1:nrow(newdata) for(d in dummies) mcd[k, d] <- mean(mcd[-k, d], na.rm = FALSE) } if(is.null(vns)){ # デフォルトは全変数 vns <- colnames(mcd)[-1] } else { # キャラクター型は列名から該当するダミー変数を特定し、置換する dc <- attr(x = r_obj$terms, which = "dataClasses")[-1] tmp <- c() for(vn in vns){ if(!exists(dc[vn]) || "numeric" == dc[vn] || "integer" == dc[vn]){ tmp <- c(tmp, vn) } else { len <- length(vn) tmp <- c(tmp, dummies[substring(dummies, 1, len) == vn]) } } vns <- tmp } # I(...)の限界効果自体は求められないので、除外する # 中の変数で評価してください vns <- excludeME(vns) # 閾値とその他の係数の持ち方を揃える if(inherits(r_obj, "polr")){ ce <- coef(r_obj) tce <- r_obj$zeta nt <- length(tce) if(interval){ replacement <- function(vcov, nov){ mt <- 1:nov nt <- (nov + 1):ncol(vcov) m <- vcov[, c(nt, mt)] m <- m[c(nt, mt), ] m } vcov <- replacement(vcov(r_obj), length(ce)) } link <- r_obj$method } else if(inherits(r_obj, "clm")) { nt <- length(coef(r_obj)) - ncol(mcd) + 1 k <- 1:nt tce <- coef(r_obj)[k] ce <- coef(r_obj)[-k] if(interval) vcov <- vcov(r_obj) link <- r_obj$link } else { stop("incorrect object!") } # 境界の係数を前に出し、説明変数の係数を負にする coef <- matrix(c(tce, -ce), 1) # 単項だけではなく、クロス項（e.g. x*z）の部分も変化させる terms = strsplit(colnames(mcd), ":") chgvalue <- function(v, data, name, value){ v2 <- numeric(length(v)) for(i in 1:length(v)){ cns <- terms[[i]] r <- 1 # I(...)の処理。交差項にI(...)は入らないようなので、ループの外に出している if(0 < length(grep("^I\\([^)]+\\)$", cns[1]))){ data[name] <- value v2[i] <- with(data, eval(parse(text = cns[1]))) next } for(j in 1:length(cns)){ # vに値が無い場合は、dataを探す r <- r*((cns[j] != name) * ifelse(is.na(v[cns[j]]), data[, cns[j]], v[cns[j]]) + (cns[j] == name) * value) } v2[i] <- r } v2 } # 限界効果を計算する変数ベクトルから、クロス項を除外しておく vns <- vns[grep("^[^:]+:[^:]+$", vns, invert = TRUE)] # vnsが空ではないことをチェック if(0>=length(vns)) stop("no effective variables in vns!") cdf <- function(a) exp(a)/(1 + exp(a)) if("probit" == link) cdf <- pnorm gpredict <- function(nd, beta = coef){ X <- cbind(diag(nt)[rep(1:nt, each = nrow(nd)), ], nd[rep(1:nrow(nd), nt), -1]) a <- beta %*% t(X) # Xβを確率密度にする関数は引数で与えた方がよいかも p <- matrix(cdf(a), nrow(nd)) r <- matrix(NA, nrow(nd), ncol(p) + 1) r[, 1] <- p[, 1] j <- 2 while(j <= ncol(p)){ r[, j] <- p[, j] - p[, j - 1] j <- j + 1 } r[, ncol(p) + 1] <- 1 - p[, ncol(p)] # 被説明変数の値ごとの確率を戻す r } lapply(1:nrow(newdata), function(i){ nd <- rbind(mcd[i, ], mcd[i, ]) t(sapply(vns, function(vn){ if(rev.dum && vn %in% dummies){ nd[1, ] <- chgvalue(nd[1, ], newdata[i, ], vn, 0) nd[2, ] <- chgvalue(nd[1, ], newdata[i, ], vn, 1) m <- gpredict(nd) me <- m[2, ] - m[1, ] } else { value <- ifelse(vn %in% colnames(nd), nd[1, vn], newdata[i, vn]) nd[1, ] <- chgvalue(nd[1, ], newdata[i, ], vn, value - h) nd[2, ] <- chgvalue(nd[1, ], newdata[i, ], vn, value + h) m <- gpredict(nd) me <- (m[2, ] - m[1, ])/(2*h) } names(me) <- levels(r_obj$model[, 1]) if(interval){ # ∂F(Xβ)/∂β を求める dFdb <- numeric(length(coef)) # nd <- matrix(mcd[i, ], 1) # ndをリセット beta1 <- beta2 <- coef # dFdbの奇数行がnd[1, ]の、偶数行がnd[2, ]をXとした∂F(Xβ)/∂β dFdb <- sapply(1:length(coef), function(j){ beta1[j] <- beta1[j] - h beta2[j] <- beta2[j] + h (gpredict(nd, beta2) - gpredict(nd, beta1))/(2*h) }) # 標準誤差を計算 # https://jslsoc.sitehost.iu.edu/stata/ci_computations/spost_deltaci.pdf の式(10) se <- sapply(1:ncol(m), function(i){ var <- t(dFdb[i, ]) %*% vcov %*% dFdb[i, ] - t(dFdb[i + 1, ]) %*% vcov %*% dFdb[i, ] - t(dFdb[i, ]) %*% vcov %*% dFdb[i + 1, ] + t(dFdb[i + 1, ]) %*% vcov %*% dFdb[i + 1, ] sqrt(var) }) a <- 1 - level margin <- qt(1 - a/2, r_obj$df) * se r <- cbind(me, me - margin, me + margin) names(r) <- c(names(me), sprintf("lwr %s", names(me)), sprintf("upr %s", names(me))) return(r) } return(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) }) model <- g ~ x1 + x2 + d library(ordinal) r_clm <- clm(model, data = df01) # summary(r_clm) # 平均値の限界効果 calcMarginalEffect(r_clm)[[1]] # 平均値の限界効果（信頼区間も計算） calcMarginalEffect(r_clm, interval = TRUE) # 平均値の限界効果（x1だけ計算） calcMarginalEffect(r_clm, vns = "x1")[[1]] # (x1, x2, d)が(0, 0, "a")と(1, 1, "b")の限界効果 calcMarginalEffect(r_clm, df01, data.frame(x1 = c(0, 1), x2 = c(0, 1), d = c("a", "b")))
クロス項や2乗項が入っていても計算します。
r_clm <- clm(g ~ x1*x2 + d, data = df01) calcMarginalEffect(r_clm)[[1]] r_clm <- clm(g ~ I(x1^2) + x1 + x2 + d, data = df01) calcMarginalEffect(r_clm)[[1]]
あっているのか間違っているのか良く分からないわけですが、ererパッケージのocMEとだいたい同じ結果になっているので、あっているということで。
library(MASS) library(erer) r_polr <- polr(model, data = df01, Hess = TRUE) ocME(r_polr, digits = 5) round(calcMarginalEffect(r_polr)[[1]], 5)
平均限界効果
戻り値が個体ごとのリスト構造になっており、平均限界効果の計算がしづらいかも知れないので、計算例を挙げて起きます。
# 平均限界効果 calcAverageMarginalEffect <- function(mes){ cns <- colnames(mes[[1]]) rns <- rownames(mes[[1]]) ame <- matrix(NA, length(rns), length(cns), dimnames = list(rns, cns)) for(i in rns){ for(j in cns){ ame[i, j] <- mean(sapply(mes, "[[", i, j)) } } ame } calcAverageMarginalEffect(calcMarginalEffect(r_clm, df01, df01))
*1:marginaleffectsは平均限界効果は簡単に計算できるのですが、平均値におけるmarginalな効果はダミー変数があるとうまく出せそうになかったです。
