Tobitモデルでの信頼区間（予測区間）

線形回帰において被説明変数の0左打ち切りを行い、0以上の値を取るようなTobitモデル適用を考えています。
右上がりの回帰直線部分の信頼区間や予測区間同様に、打ち切りのy=0部分にも（意味のある）信頼区間や予測区間を考えることは可能なのでしょうか。
打ち切り部分に同区間を記載した作図を見たことがなく。。

No.5 ベストアンサー 回答者： kamiyasiro 回答日時： 2024/02/06 16:20

No.4です。

見つけました！
作図例が無いことも、納得しました。


No.4で示したサイトのリンク先に、

『Confidence intervals for Tobit model in package AER in R』
https://stats.stackexchange.com/questions/169184 …

があります！

ここに書いてあるように、predict()関数で、予測区間も出力されるようです。

ただし、次のような注意書きがありました。まさにご質問者様が懸念されているお話ですね。

①は何を言っているかというと、y方向の誤差としては正しいが、x方向の誤差（どこで打ち切るかという値の誤差範囲）を正しく言っているのではない、ということです。
従属変数ｙが誤差を持つ式ですから当然ですね。
でも、幅を描いてしまうと、打ち切りの閾値の誤差範囲（横幅）だと読み取ってしまう人が出て来そうですね。だから作図例がないのでしょう。

～～～～～～～～～～～～～～～～～～～～～～

①「Rでの打ち切り回帰（リンクあり）」で前述したように、この手法は、打ち切り変数の期待値ではなく、潜在的な打ち切りなし変数の期待値を予測します。したがって、単に 0 でカットすることは理想的な解決策ではない可能性があります。predict() y.star y

②予測区間について:予測の推定標準誤差を取得する場合と同様に設定できます。predict(fit, se.fit = TRUE) predict.glm()

～～～～～～～～～～～～～～～～～～～～～～

私も大変勉強になりました。
特に、横幅を誤解する恐れがあることについて、教訓を得ました。
ありがとうございました。

この回答へのお礼

ご回答ありがとうございます。
なるほど、信頼区間・予測区間を記載するなら、ひとまずは回帰直線部分に限っておくことにします。

お礼日時：2024/02/08 05:16

No.4 回答者： kamiyasiro 回答日時： 2024/02/06 15:53

No.2です。

回帰線の信頼区間は定石通りブートストラップでやると提案した者です。

その後調べていたら、トービットモデルの「回帰係数の信頼区間」（回帰線の信頼区間ではなく）を求める方法↓が出ているのを見つけました。（Rですが・・・）

https://stats.stackexchange.com/questions/338360 …

だったら、回帰係数をその区間で正規乱数で発生させ、大量のモンテカルロシミュレーションを行って大量の回帰線を引けば、回帰線の信頼区間が出るかな、と考えました。手間は同じですが・・・。
（というか、回帰係数の信頼区間も、実はブートストラップだったりして）

なお、予測区間については、この信頼区間に「回帰残差から計算される偶然誤差」を分散の加法性を使って乗せることになりますが、この誤差モデルについては、いまだ調査中です。

等分散モデルで良いのでしょうか。
でも、負の予測区間が出来てしまいますよね。

ところで、ご質問者様は「y=0部分にも信頼区間や予測区間を考えることは可能か」とご質問されていますが、それ以外の箇所は信頼区間や予測区間を求めた、ということですか？

この回答へのお礼

ご回答ありがとうございます。
> ところで、ご質問者様は「y=0部分にも信頼区間や予測区間を考えることは可能か」とご質問されていますが、それ以外の箇所は信頼区間や予測区間を求めた、ということですか？
最尤法では未算出です。打ち切り考慮なしでOLSでは算出しました。等分散を仮定していますので、そうですね、打ち切っても予測区間が負数を含みますよね。。

お礼日時：2024/02/08 05:12

No.3 回答者： stomachman 回答日時： 2024/02/06 12:21

> 降水量などのイメージかな

そういうのには、magnitudeを使いますね。降水量の対数を取ったのをデータとする。だから、降水量0というのはそもそもデータのうちに数えない。たとえば「千年に一度の大豪雨」なんて言っているのは、「降水量rの対数と、時間あたりの頻度の対数との両対数グラフを描くと大体直線になるよ」という経験則に基づいているわけで、これは
　　log(r) = a log(t) + b + ε
というモデル。だから
　　r = (t^a) exp(b + ε)
を考えていることになる。

　ちなみに、これと似ていても、「基準値Cからの増分rの対数と、説明変数tの対数をプロットすると大体直線になる。ただし、残差εはtやrとは無関係な外乱によって生じている」ということだと、
　　r = (t^a) exp(b) + ε
というモデルであり、rは0以下になることもあるんで、正直に非線形モデルとして扱う。こちらのモデルは残差を「絶対誤差」として扱っているのに対し、大豪雨のモデルは残差を「rに対する相対誤差」として扱っているんで、rが大きいほど大きな残差を許している。このように、残差の発生メカニズムも考慮しないとモデルが決まらない。（もちろん、「信頼区間」を心配する以上は、当然そこも考えてあるはずですけれども。）

　ところで、非線形モデルのパラメータ推定は反復計算をやるんで、その出発値を決めるのに線形モデルによる近似を考えるのは意味がある。すなわち、rが0以下であるようなデータを全部無視した上で、
　　a log(t) + b = log( r - ε) = log(r) + log( 1 - ε/r) ≒ log(r) - ε/r
より、
　　r log(r) = a log(t) r + b r + ε
という線形回帰をやれば出発値aが簡単に決められる。というわけで、非線形モデルにおいても、線形回帰は使いどころがあります。で、「rが0以下であるようなデータを全部無視」というところがTobitっぽいけれども、飽くまでも反復計算の出発値の決定をやっているに過ぎないんだから、この段階で残差の分布を心配しても意味がない。

> 被説明変数は「プログラムのバグ数」ベースの値

　「小中学校の授業で児童100人を対象に、全員に同じお題を与えてプログラム作らせた。それぞれの児童について、いくつバグがあったか」という話ですら、バグの数を数えるのは困難でしょう。プログラムの体をなしていないシッチャカメッチャカな作品の「バグの数」なんて定義できないからです。やれるとすれば「各児童の作品が、予め用意したテストをパスするかどうか」の判定でしょう。すると、ロジスティックモデル（被説明変数はテストにパスする確率、説明変数はたとえば、国語の成績とか）を使うしかなかろう。
　一方、実用的な大規模なプログラムのバグ数を数えるのも無理。テストして探してたんでは「バグはこれで全部」と宣言することは不可能で、program verification（プログラムの仕様を論理式で書いて、これをプログラムが満たすことを(数学の意味で）証明する）によるしかないのだが、実際にこれをやったという話は寡聞にして知らない。大規模なプログラムでは、そもそも「仕様を論理式で書」くのが大変だからです。
　すると、もしかして「プログラムのバグ数」ではなく、「あるプログラムについて、（これまでに）見つかったバグ数」の話をなさっているのではないか。
　そうだとすると、さて、線形回帰になっていて、かつナンセンスでないようなモデルなんてものが果たしてあるんだろうか。
　たとえばテストに掛けたman・power xについて、
　　N = B(1 - exp(-x/A))
のようなモデルはよく使われる。（直接観測できない）「プログラムのバグ数」がパラメータBで、Aは「テストの効率の悪さ」的なものを表す。実測値N(x)がこのモデルによくfitする場合、パラメータを推定すると「いくつぐらいバグが残っているか」が推測できるからです。しかしこれは非線形のモデルだなあ。

…てなことを思ってしまう。どうもご質問の情報だけでは、寄り添って考えるのは難しいが、ま、いずれにしてもTobitとはあんまり関係なさそうな気はする。

この回答へのお礼

ご回答ありがとうございます。
> 「プログラムのバグ数」ではなく、「あるプログラムについて、（これまでに）見つかったバグ数」
おっしゃる通り後者を構想していました。
手元のデータでは、一旦は線形関係が見出せていそうなのであとは0近傍の処理かと思っていましたが、、既存モデルの調べが足りていなかったので確認してみることにします。
詳細な回答を頂きまして、まわりに有識者が少ないこともあり、知見が広がり助かります。

お礼日時：2024/02/08 05:12

No.2 回答者： kamiyasiro 回答日時： 2024/02/04 23:19

ロジスティック回帰でも、ｙ＝０や１の部分に信頼区間の幅がありますので、トービットモデルでも回帰線がｙ＝０になる時点でプラス側に信頼区間や予測区間がありますよね。

私はＲ使いなんですが、
ＲではAERパッケージでトービット回帰はできるようですが、predict()関数でconfidence intervalsが算出できるのかどうかは知りません。

たとえできなくても、OLS以外の信頼区間算出はブートストラップ法でできるので、私だったらブートストラップでやります。
プログラムを書かなければなりませんが・・・。

予測区間については方策が浮かびません。
誤差εが等分散仮定なら２乗和で出来そうですが、ｙ＝０でε＝０に収束するのですよね。それに飽和状態もε＝０ですよね。

誤差モデルはbinomialなんでしょうか？

この回答へのお礼

ご回答ありがとうございます。
誤差はbinomial想定です。

お礼日時：2024/02/08 05:12

No.1 回答者： stomachman 回答日時： 2024/02/04 20:58

Tobitと仰るからには「被説明変数が負の値の場合は観測できない」ってことかしらん。

その場合、観測できなかった（負に違いない）被説明変数の値がどんな分布に従うかがわからなければ、もちろん区間推定などできない。しかし、もし分布のパラメータの値まで含めて理論的に分かっているのなら、当然、それを使って区間を推定できる（そんな状況は滅多にないだろうけど）。しかし、被説明変数が非負である部分についてはフツーの推定は観測したデータの経験分布（正規分布などで近似して扱うかどうかはどっちでもよくて、ともかくそれ）を使って行うんだとしますと、両者のやり方が全然違うので、少なくともこれらをごっちゃにしてプロットするのはおかしいでしょう。
　一方、もし「打ち切り」と仰るのが、「被説明変数が負の値の場合も観測できているにもかかわらず、回帰モデルのパラメータを推定するにあたって、被説明変数が負の値を全部"データの欠落"とみなして強制的にTobitで扱う」という話だとしますと、こんな変なことはまずやらない（∵ 推定した回帰モデルにバイアスが入る）んだけれど、「いや、負の時には観測がメチャクチャ不安定でデータがまるで信用できんのだ」というような事情があれば、ま、分からなくもない。被説明変数が負の値の部分についてもデータがあるんですから、経験分布を使ってフツーに区間推定をすれば答は得られる。しかし、もしそのデータがまるで信用できんのであれば、それを使って推定した区間なんて、わざわざやるだけの意味はない。

この回答へのお礼

詳細なご回答ありがとうございます。理解が進みました。
被説明変数は「プログラムのバグ数」ベースの値（非負の実数）でして、「負の値の場合は観測できない」と言うよりは、説明変数が一定以下では、多くは0（しばしば正の実数）が観測されている、という状況になります。
降水量などのイメージかなと思います。0も実際に観測される値と言う意味ではtwo-part modelなどがより適するのかもしれませんが、一旦Tobitで考えていたところでした。

お礼日時：2024/02/06 07:30