タイプ1トービットモデル
UPDATE: 2023-04-22 14:46:46
はじめに
ここではタイプ1トービットモデルについて、モデルの内容、モデルの推定法(最尤法)、Rでの実行方法をまとめている。タイプ2トービットモデルについてはここでは扱っていないが、下記の参考にした書籍には記載があるので、そちらを参照のこと。
トービットモデル
トービットモデルは、データの上限、下限を持っている打ち切りデータ(CensoredData)に利用される手法の1つ。目的変数が打ち切りデータの場合、通常の回帰モデルでパラメタを推定しても、最小二乗法の推定量は普遍性も一致性ももってない。
打ち切りデータの例として、小売の会員の購買金額は0よりも小さくなることはないし、店舗への会員訪問回数も同じである。なぜ打ち切りデータに通常の回帰モデルが不適なのか、直感的な説明をすると、データ全体に対してモデルの回帰係数を推定した結果、パラメタが正なり負なり推定されるわけだが、売上が0円の会員(ここでは登録したが買ってないケース)もいるわけで、これらはパラメタや説明変数がどうであれ、目的変数は0で一定であるためである。
この問題を回避するために売上が0ではないデータに限定して、モデルを解釈する方法もあるように見える。顧客の属性や行動履歴を用いて、モデリングした結果、当初の分析目的である会員全体を対象にしている分析でもなく、特定の集団を分析しているため、このモデルの解釈は限定的になる。その点を無視して、会員全体に分析結果を反映させるとセレクションバイアスが発生していることになる。売上0より大きい会員は特定の特徴を持っている可能性があるためである。
他にも、この記事を書いている時期が野球の世界大会であるWBCとかぶっているので、野球の例であれば、出場していない選手の空振り数やヒット数などのパフォーマンスは、出場していないのであれば0である。このようなケースで出場している選手のデータに限定すると、セレクションバイアスが発生する。ほかにも、一定の世帯年収がなければ教育費は0になっているデータに対して、教育費が0を除外したデータだけで分析するわけにもいかない。
モデルの推定法
ここでは、個人を\(i\)、個人が持つ説明変数を\(x_{i} = (x_{i1},x_{i2},...,x_{ik})^{t}\)とする。トービットモデルでは、目的変数\(y_{i}\)に対して、\(y_{i}\)に対する効用を表す潜在変数\(y^{*}_{i}\)を定義する。\(y^{*}_{i}\)は観測されないが、観測できる\(y_{i}\)とは下記のように関係を考える。
\[ \begin{eqnarray} y^{*}_{i} = \begin{cases} y^{*}_{i} & ( y^{*}_{i} \gt 0 ) \\ 0 & ( y^{*}_{i} \le 0 ) \end{cases} \end{eqnarray} \]
$y^{*}_{i} \(では\)y^{*}_{i}\(と\)y^{*}_{i}\(が一致し、\)y^{*}_{i} \(では\)y^{*}_{i}\(は0または負の値を持つことになる。効用を表す潜在変数\)y^{*}_{i}\(を利用することで、説明変数が大きくなると\)y^{*}_{i}\(が大きくなる、説明変数が小さくなると\)y^{*}_{i}$は小さくなる、という自然な仮定を考えることができる。これであれば通常の回帰モデルを利用することは問題にならない。
\[ y^{*}_{i} = x^{t}_{i}\beta + u_{i}, u_{i} \sim N(0, \sigma^{2}) \]
ただ、\(y^{*}_{i}\)は観測できないので、\(y_{i}\)を使って分析すると、パラメタの普偏性や一致性が損なわれ、データを限定しても同様である。そのため、\(y^{*}_{i}\)を使ってモデルを計算できるように工夫してパラメタを推定する必要がある。トービットモデルの推定方法としては最尤法とヘックマン2ステップ推定法があるらしいが、後者は勉強不足でわからないので、最尤法を利用する。
まず、領域を分けて考えると\(y_{i} \gt 0\)の領域では、\(y^{*}_{i}\)と\(y^{*}_{i}\)が一致するので、\(y_{i}\)を観測する確率は下記のように書ける。ここで\(\phi\)は標準正規分布を表す。
\[ \begin{eqnarray} P(y_{i}) &=& P \left(y^{*}_{i} \right) \\ &=& P \left(x^{t}_{i}\beta + u_{i} \right) \\ &=& P \left( x^{t}_{i}\beta + \epsilon_{i} \sigma \right)\\ &=& \frac{1}{\sigma} \phi \left( \frac{ y^{*}_{i} - x^{t}_{i}\beta}{\sigma} \right) \\ \end{eqnarray} \]
次に、\(y_{i} = 0\)のとき、\(y_{i}\)を観測する確率を考える。さきほどの\(u_{i}\)に対して、\(\epsilon_{i} = \frac{u_{i}}{\sigma}\)とすると、\(\epsilon_{i}\sim N(0, 1)\)となり、標準正規分布に従うことがわかる。\(\Phi\)は累積分布関数。
\[ \begin{eqnarray} P(y_{i}=0) &=& P \left(y^{*}_{i} \le 0 \right) \\ &=& P \left(x^{t}_{i}\beta + u_{i} \le 0 \right) \\ &=& P \left( \frac{u_{i}}{\sigma} \le - \frac{x^{t}_{i}\beta}{\sigma} \right)\\ &=& P \left( \epsilon_{i} \le - \frac{x^{t}_{i}\beta}{\sigma} \right)\\ &=& 1 - \Phi\left( \frac{x^{t}_{i}\beta}{\sigma}\right)\\ \end{eqnarray} \]
これで各領域で確率を計算できそうなので、組み合わせると尤度関数が構築できる。\(\delta(y_{i})\)は指示関数で、領域によって振り分けを行う。
\[ L(\beta, \sigma^{2}|{y_{i},x_{i}}) = \displaystyle \prod_{i=1}^{n} \left[\frac{1}{\sigma} \phi \left( \frac{ y_{i} - x^{t}_{i}\beta}{\sigma} \right)\right]^{1-\delta(y_{i})}\left[ 1 - \Phi\left( \frac{x^{t}_{i}\beta}{\sigma}\right)\right]^{\delta(y_{i})} \]
参考書を読んでいるときに使った手書きのメモはこれ。デジタルにする気力がなかったので、そのままアップしておく。
Rで最尤法を使って計算してみる。まずはパラメタが正しく推定できているか、ダミデータを作成する。\(\alpha=1, \beta=2\)とする。また、\(y \lt 0\)のとき、打ち切りデータとなり、0になる。この打ち切りデータに対して最小二乗法の回帰モデルで推定すると、うまく推定できていないことがわかる。
set.seed(1989)
N <- 5000
x <- rnorm(N)
e <- rnorm(N)
a <- 1
b <- 2
y <- a + b * x + e
y_cens <- ifelse(y < 0, 0, y)
df_cens <- data.frame(x, y, y_cens)
fit <- lm(y_cens ~ x, data = df_cens)
fit##
## Call:
## lm(formula = y_cens ~ x, data = df_cens)
##
## Coefficients:
## (Intercept) x
## 1.482 1.359データを限定しても、正しく推定できていない。
fit_filter <- lm(y_cens ~ x, data = dplyr::filter(df_cens, y>0))
fit_filter##
## Call:
## lm(formula = y_cens ~ x, data = dplyr::filter(df_cens, y > 0))
##
## Coefficients:
## (Intercept) x
## 1.415 1.653これに対してトービットモデルでパラメタを推定してみる。ここでは対数尤度関数に変換して利用している。推定された結果を見る限り問題なく推定されているように見える。
# 初期値
p <- c(sigma = 1, a = 1, b = 1)
X <- model.matrix(y_cens ~ x, data=df_cens[,c(1,3)])
y <- df_cens$y_cens
optim(par = p,
fn = function(p) {
sigma <- p[1]
beta <- matrix(p[-1], ncol = 1)
xb <- X %*% beta
-1 * sum(ifelse(y <= 0, log(1 - pnorm(xb/sigma)), log(dnorm((y - xb)/sigma)/sigma)))
},
method = "BFGS",
hessian = TRUE)$par## sigma a b
## 0.998270 1.001899 2.019256Rでの実践
参考にしている書籍のデータをお借りして、トービットモデルへの理解を深める。タイプ2トービットモデルはここでは扱わないが、書籍では説明されているので、そちらを参照。
データセットの内容は下記の通り。オフラインの店舗の企業が、オンライン店舗を出店し、顧客分析を行う際の例。会員全体を対象にしている分析で、オンライン店舗での購買金額が、顧客の属性や行動履歴によって、どのように説明されているかを知りたいとする。それがわかれば、似たような会員をオンライン店舗での購入につなげる施策を検討できる。打ち切りデータなので、オンラインで購入していない会員のオンライン店舗の売上は0となっている。
| カラム名 | 内容 |
|---|---|
CID |
会員ID |
Purchase |
会員のオンライン店舗の利用有無(1=利用あり) |
Amount |
会員のオンライン店舗の購買金額 |
Sex |
会員の性別(1=男性) |
Age |
会員の年齢 |
Fsize |
会員の家族人数 |
Income |
会員の対数収入 |
Ownhouse |
会員の持ち家有無(1=あり) |
Crossbuying |
会員の購入カテゴリ数 |
Pfreq |
会員の購買頻度 |
Rduration |
会員の入会期間(月数) |
トービットモデルはcensRegパッケージのcensReg()関数で利用できる。IncomeとOwnhouse以外は有意となっている。このモデルでは、男性で家族の人数が多いほうがオンライン店舗での購入金額は大きくなり、年齢では若い方が購入金額が大きくなることがわかる。
library(censReg)
data_chap7 <- read.csv('~/Desktop/chapter_07.csv')
fit_tobit <-
censReg(Amount ~ Sex + Age + Fsize + Income + Ownhouse + Crossbuying,
data = data_chap7)
summary(fit_tobit)##
## Call:
## censReg(formula = Amount ~ Sex + Age + Fsize + Income + Ownhouse +
## Crossbuying, data = data_chap7)
##
## Observations:
## Total Left-censored Uncensored Right-censored
## 1000 689 311 0
##
## Coefficients:
## Estimate Std. error t value Pr(> t)
## (Intercept) -21.01248 24.85756 -0.845 0.398
## Sex 8.88941 0.97833 9.086 < 2e-16 ***
## Age -0.46155 0.03943 -11.705 < 2e-16 ***
## Fsize 4.20824 0.56364 7.466 8.26e-14 ***
## Income 0.11736 1.63300 0.072 0.943
## Ownhouse 0.99016 1.14626 0.864 0.388
## Crossbuying 1.15683 0.08207 14.096 < 2e-16 ***
## logSigma 2.36784 0.04394 53.887 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Newton-Raphson maximisation, 7 iterations
## Return code 1: gradient close to zero (gradtol)
## Log-likelihood: -1377.252 on 8 Dfトービットモデルのパラメタは解釈に注意が必要で、通常の最小二乗法の回帰モデルは、説明変数の\(E(y_{i})\)に対する限界効果である一方で、トービットモデルは\(E(y^{*}_{i})\)に対する説明変数の限界効果である。もし最小二乗法の回帰モデルの係数と比較するのであれば、トービットモデルでの\(E(y_{i})\)に対する説明変数の限界効果を知る必要がある。これは下記の式で計算できる。
\[
\frac{\partial E(y_{i})}{\partial x_{ji}} = \beta_{j} \Phi \left(
\frac{x^{t}_{i} \beta}{\sigma} \right)
\] この場合、限界効果は\(x_{i}\)に依存するため、平均を使うことで、平均限界効果を計算する。これはmargEff()関数で計算できる。
margEff(fit_tobit)## Sex Age Fsize Income Ownhouse Crossbuying
## 1.89474002 -0.09837827 0.89696761 0.02501544 0.21104787 0.24657241平均限界効果を用いると、最小二乗法の回帰モデルと係数の比較ができる。平均はmargEff()の中で、censRegクラスのオブジェクトがすでにもっている。
# purrr::map_dbl(data_chap7,mean)でも同じ
fit_tobit$xMean## (Intercept) Sex Age Fsize Income Ownhouse
## 1.00000 0.34700 43.50000 2.93400 15.16593 0.37000
## Crossbuying
## 12.98400この関数の中身は下記の通り。
getS3method("margEff", "censReg")## function (object, xValues = NULL, vcov = NULL, calcVCov = TRUE,
## returnJacobian = FALSE, vcovLogSigma = TRUE, ...)
## {
## allPar <- coef(object, logSigma = FALSE)
## isPanel <- "sigmaMu" %in% names(allPar)
## if (isPanel) {
## stop("the margEff() method for objects of class 'censReg'",
## " can not yet be used for panel data models")
## }
## sigma <- allPar["sigma"]
## beta <- allPar[!names(allPar) %in% c("sigma")]
## if (is.null(xValues)) {
## xValues <- object$xMean
## if (length(xValues) != length(beta)) {
## print(beta)
## print(xValues)
## print(object$xMean)
## stop("cannot calculate marginal effects due to an internal error:",
## " please contact the maintainer of this package")
## }
## }
## else if (!is.vector(xValues)) {
## stop("argument 'xValues' must be a vector")
## }
## else if (length(xValues) != length(beta)) {
## stop("argument 'xValues' must be a vector with the number of elements",
## " equal to the number of estimated coefficients without sigma (",
## length(beta), ")")
## }
## else {
## names(xValues) <- names(object$xMean)
## }
## xBeta <- drop(crossprod(xValues, beta))
## zRight <- (object$right - xBeta)/sigma
## zLeft <- (object$left - xBeta)/sigma
## result <- beta[!names(beta) %in% c("(Intercept)")] * (pnorm(zRight) -
## pnorm(zLeft))
## names(result) <- names(beta)[!names(beta) %in% c("(Intercept)")]
## if (calcVCov || returnJacobian) {
## jac <- matrix(0, nrow = length(result), ncol = length(allPar))
## rownames(jac) <- names(result)
## colnames(jac) <- names(allPar)
## for (j in names(result)) {
## for (k in names(allPar)[-length(allPar)]) {
## jac[j, k] <- (j == k) * (pnorm(zRight) - pnorm(zLeft)) -
## (beta[j] * xValues[k]/sigma) * (dnorm(zRight) -
## dnorm(zLeft))
## }
## jac[j, "sigma"] <- 0
## if (is.finite(object$right)) {
## jac[j, "sigma"] <- jac[j, "sigma"] - (beta[j]/sigma) *
## dnorm(zRight) * zRight
## }
## if (is.finite(object$left)) {
## jac[j, "sigma"] <- jac[j, "sigma"] + (beta[j]/sigma) *
## dnorm(zLeft) * zLeft
## }
## }
## if (calcVCov) {
## if (is.null(vcov)) {
## vcov <- vcov(object, logSigma = FALSE)
## }
## else {
## errMsg <- paste0("argument 'vcov' must be a symmetric ",
## ncol(jac), " x ", ncol(jac), " matrix")
## if (!is.matrix(vcov)) {
## stop(errMsg)
## }
## else if (!isSymmetric(vcov)) {
## stop(errMsg)
## }
## else if (nrow(vcov) != ncol(jac) || nrow(vcov) !=
## ncol(jac)) {
## stop(errMsg)
## }
## if (vcovLogSigma) {
## if ("sigma" %in% c(rownames(vcov), colnames(vcov))) {
## warning("Please make sure that argument 'vcov' includes",
## " the covariances of 'log(sigma)' rather than those of 'sigma'",
## " or set argument 'vcovLogSigma' to 'FALSE'")
## }
## vcov[nrow(vcov), ] <- vcov[nrow(vcov), ] *
## allPar["sigma"]
## vcov[, nrow(vcov)] <- vcov[, nrow(vcov)] *
## allPar["sigma"]
## }
## else {
## if ("logSigma" %in% c(rownames(vcov), colnames(vcov))) {
## warning("Please make sure that argument 'vcov' includes",
## " the covariances of 'sigma' rather than those of 'log(sigma)'",
## " or set argument 'vcovLogSigma' to 'TRUE'")
## }
## }
## }
## attr(result, "vcov") <- jac %*% vcov %*% t(jac)
## }
## if (returnJacobian) {
## attr(result, "jacobian") <- jac
## }
## }
## attr(result, "df.residual") <- object$df.residual
## class(result) <- c("margEff.censReg", class(result))
## return(result)
## }
## <bytecode: 0x10a2e3c40>
## <environment: namespace:censReg>