信頼区間とはなんぞやというのをメモしておこうと思って、簡単なデータで回帰分析を行って図をつくろうかと思ったら、Rのpredict()関数の使い方に落とし穴があったので復習がてらメモ……。
とりあえず単回帰分析する
Rの練習用データセット「cars」をつかいます。*1
車のスピードと制動距離（or 停止距離）ですかね。
> head(cars) # Rの練習用データセット「cars」の中身 speed dist 1 4 2 2 4 10 3 7 4 4 7 22 5 8 16 6 9 10
相関係数と散布図をみておきます。
> cor(cars$speed, cars$dist) [1] 0.8068949 > plot(cars)
とりあえず回帰分析します。
# lm()で回帰分析 > cars.lm <- lm(dist ~ speed, data=cars) > summary(cars.lm) Call: lm(formula = dist ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -29.069 -9.525 -2.272 9.215 43.201 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.5791 6.7584 -2.601 0.0123 * speed 3.9324 0.4155 9.464 0.00000000000149 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.38 on 48 degrees of freedom Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438 F-statistic: 89.57 on 1 and 48 DF, p-value: 0.00000000000149
散布図に回帰直線をひきます。
> # 散布図に回帰直線を引く > plot(cars) > abline(cars.lm, + lwd=2, + col="red" + )
べつに意味は無いですが、predict()関数で予測値を計算してその折れ線グラフを書けば、当然、回帰直線と同じ線が引けます。x軸の範囲が限定されますが。
> # predict()関数で予測値を得る > cars.predict <- predict(cars.lm) > > # 予測値の中身をみるとベクトルになっている。 > head(cars.predict) # head()は頭から6要素分（6行分）のデータを確認する関数 1 2 3 4 5 6 -1.849460 -1.849460 9.947766 9.947766 13.880175 17.812584 > > # 元のデータのプロット > plot(cars, + xlim=c(min(cars$speed)-1, max(cars$speed)+1), # 軸の目盛をムダに丁寧に調節ｗ + ylim=c(min(cars$dist)-1, max(cars$dist)+1) + ) > > # 予測値を折れ線でプロット > par(new=TRUE) > plot(cbind(cars[1], cars.predict), + xlim=c(min(cars$speed)-1, max(cars$speed)+1), + ylim=c(min(cars$dist)-1, max(cars$dist)+1), + axes=F, ylab="", xlab="", # par(new=T)で重ねるときは後のグラフのメモリとラベルを非表示に + type="l", + lwd=2, + col="red" + ) >
信頼区間をみてみる
次に、信頼区間を求めます。
> # 信頼区間をみるにはpredict関数を用い、intervalを"confidence"に設定する > > cars.conf <- predict(cars.lm, interval="confidence", level=0.95) > cars.conf <- as.data.frame(cars.conf) # データフレームにしとく > > > # 中身をみてみる > head(cars.conf) fit lwr upr 1 -1.849460 -12.329543 8.630624 2 -1.849460 -12.329543 8.630624 3 9.947766 1.678977 18.216556 4 9.947766 1.678977 18.216556 5 13.880175 6.307527 21.452823 6 17.812584 10.905120 24.720047 >
この「fit」ってのは回帰直線を表すデータで、「lwr」は信頼区間の下側の限界を表すデータ、「upr」は上側です。predict()の引数level=0.95に設定しているので、「95%信頼区間」を求めています。
これをさっきのグラフに重ねていきます。
# もとのデータの散布図 plot(cars, xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1) ) # 回帰直線に該当する予測値の折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.conf$fit), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1.5, lty=1, # 実線の場合は「1」 col="red" ) # 信頼区間の上側の限界値を表す折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.conf$upr), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1, lty=1, # 実線の場合は「1」 col="blue" ) # 信頼区間の上側の限界値を表す折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.conf$lwr), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1, lty=1, # 実線の場合は「1」 col="blue" )
赤い線が回帰直線で、青い線に挟まれた範囲が信頼区間です。
上のコードでは、線１本１本についていちいち表示を調節するための引数を与えているのですが、あとでメモするようにmatplot()を使えば複数の線の設定をいっぺんにできるので、上のようにごちゃごちゃ書く必要はないですね。
で、統計の勉強の初歩レベルでは定番のネタだと思うのですが、「95%信頼区間」というのは、「この範囲に真の値（母数）がある確率が95%」ということを表しているのではないという点に、注意が必要なわけですね。
信頼区間の解説については↓このページを見ればいいと思うのですが、
ここで言われてるのは要するに、95%信頼区間とは、その上側の限界が真のパラメータであった場合も、下側の限界が真のパラメータであった場合も、そのパラメータを持つモデルから得られるデータのうち尤もらしい方から95%の範囲内に、今回の計測で得られているサンプルデータ（から推定されるパラメータ）が含まれるような区間ということですね。
文章で表現するとめっちゃ回りくどいですが。
他にもたとえば「真の平均値」について考えているのであれば、「この範囲に『真の平均値』（母数）があると仮定すれば、その範囲内のどの値を『真の平均値』とするモデルであっても、そのモデルから得られるサンプルの平均値の95%区間（そのモデルの誤差分布上で、確率の低いほうから積分して5%分になる範囲を除いた区間）に、今回得られているサンプルの平均値が含まれるようになるような、『真の平均値』の範囲」という意味になると思います。
上のグラフでいうと、この範囲に収まる直線であれば、その直線が真のモデルであったときの誤差分布上の95%区間に、今回の回帰直線が含まれるということですね。たぶん。
しかしこういう言い方だけだと、どちらかというと信頼区間の「決め方」しか説明していなくて、その「本質的な意味」の説明にはなっていないですね。
これについてYahoo!知恵袋で質問しました。
[統計]「信頼区間」の概念をどう理解すれば良いかについての質問です。 - ... - Yahoo!知恵袋
（回答だけでなく返信でいろいろ続いていってます。）
回答者の方に丁寧に解説して頂いた結果、頭がだいぶ整理できました。（やはり「詳しい人に訊く」のは勉強する上で大事だ。）
私は何が知りたかったのかというと、95%信頼区間を求めると「何がうれしいのか」ということです。もっと具体的に言えば、「何が95%だから嬉しいのか」ということです。これは、さっきのページの解説だけだと、よくわからない*2。奥村先生の解説は「よくある誤解」をただすことを主目的に書かれてるんでしょう。
回答者の方も言われるように、「サンプリングして信頼区間を出す、というプロセスを１００回繰り返すと、大体９５回分の信頼区間は真の値を含む」ということなんですよね。
回答からリンクが貼られている解説サイトがあって、
5. 正規母集団からの標本に基づく推論
ここに出てくる図がわかりやすいのですが、要するに教科書に載ってるような計算方法（ルール）で信頼区間を定義して、同じ母集団から100回ぐらいサンプリングを繰り返しつつ毎回その区間を求めてみると、95回ぐらいは、その区間のなかに「真の値」が含まれるのであると。もちろん、誤差分布に関する仮定（たとえば「正規分布で、分散は○○である」みたいな）が正しければの話ですが。
ある１回のサンプリングから求めた推定値を中心に広がる「今回得られた95%信頼区間」の中に“何かがある”とイメージしてしまうと、誤解のもとになるんでしょう（ベイズ統計学の信用区間の場合はそういうイメージなんだけど）。そうではなく、「そのルールで区間を求めることにしておきさえすれば、何回もサンプリングと区間推定を繰り返した場合に、95/100ぐらいは当たる」のだということです。あくまで母数というのは１つの決まった値しかあり得ないはずで、「この区間に母数が含まれる確率」を議論することは（頻度主義的には）意味がなく、サンプリング誤差によって右に行ったり左に行ったり（上に行ったり下に行ったり）するのは、「計算された信頼区間」の方なので。
まぁ考えようによっては同じことを言ってる気もしますが、一応、Wikipediaの解説とかを見ても、区別されることになってるらしい。
言い換えると、「今回のサンプルデータから得られたこの信頼区間」のなかに母数（真の値）が含まれる確率が95%なのではなく、この「信頼区間算出ルール」の「勝率」が95%になるのだと理解しないといけないわけですね。「算出ルールの勝率」という概念で、割とうまく頭が整理されるのではないかと思います。
「信頼区間」というのは「データを入れると信頼区間が出てくるルール」なんです．アルゴリズムとか関数といってもいい．推定値ってのはみんなそうじゃないか，っていわれそうだけど，信頼区間については特にそれを強調する価値があると思う．
こういう説明もあったのですが、まさにこの「ルール」の「勝率」が「95%」であるという理解の仕方が大事だと思います。
ただ、知恵袋の回答者の方も言われるように、ベイズ主義のほうが思想がシンプルで良い気はする……。頻度主義は直観的には理解しやすくない。
……で、本来、この「信頼区間とは何か」ってのをメモっておこうと思っただけだったのですが、ついでに「予測区間」も求めておこうと思った時に、predict()関数の使い方をミスっていたので、以下メモしておきます。
予測区間も求めてみる
信頼区間に加えて予測区間も求めて描画してみます。信頼区間というのはあくまで推定されるモデルのパラメータの範囲について言っているもので、その範囲にあるモデルから得られるデータというのは、さらにサンプリング誤差を伴ってくるので、信頼区間よりも予測区間のほうが広くなります。
> cars.pred.int <- predict(object=cars.lm, # objectって普通は書かないけど後の説明上わかりやすいように + newdata=data.frame(speed=cars$speed), # ★ココに注意！ + interval="prediction", + level=0.95) > cars.pred.int <- as.data.frame(cars.pred.int) # データフレームにしとく > head(cars.pred.int) fit lwr upr 1 -1.849460 -34.49984 30.80092 2 -1.849460 -34.49984 30.80092 3 9.947766 -22.06142 41.95696 4 9.947766 -22.06142 41.95696 5 13.880175 -17.95629 45.71664 6 17.812584 -13.87225 49.49741 >
ここで注意点があって、predict()関数の第2引数の「newdata」ってのは説明変数を与える部分（説明変数がどういう値だったときの予測値を知りたいかを出力するわけだから）なんですが、ここには、
①データフレーム型であり、かつ、
②objectに与えたモデルの説明変数の名前を列名称とする列が存在する
ようなデータを与えないとエラーが出てしまいます。
なんか、いろいろ適当に書いてみて動かなかったので、「あれ？何か間違ってたっけ？」と思って調べたら、たとえば、
lm - R: numeric 'envir' arg not of length one in predict() - Stack Overflow
このページで説明されてました。
ちなみに、予測区間じゃなくて信頼区間を求めた時（intereval="confidence"の時）は、そもそも第2引数newdataには何も入れなくても結果が出ました。しかし予測区間を求めるときは与えないとダメらしい。
ちなみに他のサイト（たとえば↓）を見ると、信頼区間を求めるときも、たとえば「new」とかいう名前のデータフレームを予めつくって、そのデータフレームに、説明変数と同じ名称をもつ列を設けていますね。
なので、必ずnewdata=を与えるようにと覚えておいたほうが良いとは思います。「newdata=」という表記自体は省略する人が多いと思いますが。
70. 単回帰分析
おっと危ない：信頼区間と予測区間を混同しちゃダメ - Take a Risk：林岳彦の研究メモ
また、私はnewdataのところに、元のデータの説明変数を与えていますが、↑のサイトのように等間隔なデータを与えてそれをX軸としたほうが良いとも思います。元データの説明変数がけっこう飛び飛びの値だったりするとグラフが汚くなりますし。そもそも、元の分析データの説明変数（x軸）はたまたま得られた数値であって、説明変数が他の値だった場合にどうなるかを考えるのが予測であるとも言えますしね。
なお、newdataに与えるデータフレームは、無関係な列を含んでいても問題ないようです。objectに入れたモデルの説明変数と同じ名称の列を探してくるようになっているようですね。↓のようなことをしても大丈夫でした。
> test <- predict(object=cars.lm, + newdata=data.frame(speed=cars$speed, muda1=seq(1, 50), muda2=rep("muda", 50)), + interval="prediction", + level=0.95) > head(test) fit lwr upr 1 -1.849460 -34.49984 30.80092 2 -1.849460 -34.49984 30.80092 3 9.947766 -22.06142 41.95696 4 9.947766 -22.06142 41.95696 5 13.880175 -17.95629 45.71664 6 17.812584 -13.87225 49.49741
さてこの予測区間を、さっきの信頼区間つきのグラフに重ねてみます。
# もとのデータの散布図 plot(cars, xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1) ) # 回帰直線に該当する予測値の折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.conf$fit), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1.5, lty=1, # 実線の場合は「1」 col="red" ) # 信頼区間の上側の限界値を表す折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.conf$upr), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1, lty=1, # 実線の場合は「1」 col="blue" ) # 予測区間の上側の限界値を表す折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.conf$lwr), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1, lty=1, # 実線の場合は「1」 col="blue" ) # 信頼区間の上側の限界値を表す折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.pred.int$upr), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1, lty=2, # 実線の場合は2 col="blue" ) # 信頼区間の上側の限界値を表す折れ線 par(new=TRUE) plot(cbind(cars$speed, cars.pred.int$lwr), xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=1, lty=2, # 点線の場合は2 col="blue" )
青い点線になっているのが、予測区間の上下の限界線です。めっちゃ広いですね。
matplot()でグラフを描いたほうが早い
ところで、上記のようにグラフを何本も書くのであれば、plot()じゃなくてmatplot()というのを使ったほうがいいですね。引数を行列として与えて、複数のグラフを一気に描くやつです。
# もとの散布図 plot(cars, xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1) ) # 直線５本をまとめる data <- cbind( cars.conf$fit, # 回帰直線 cars.conf$upr, # 信頼区間の上側 cars.conf$iwr, # 信頼区間の下側 cars.pred.int$upr, # 予測区間の上側 cars.pred.int$lwr # 予測区間の下側 ) par(new=TRUE) matplot( x = cars$speed, # 説明変数 y = data, # まとめたデータ xlim=c(min(cars$speed)-1, max(cars$speed)+1), ylim=c(min(cars$dist)-1, max(cars$dist)+1), axes=F, ylab="", xlab="", type="l", lwd=c(1.5, 1, 1, 1, 1), lty=c(1, 1, 1, 2, 2), col=c("red", "blue", "blue", "blue", "blue") )
同じ図がかけました。
