データと回帰曲線 †回帰分析で推定された曲線を回帰曲線といいます。単回帰出力結果の読み方の例でいうと、 y = -6.0930 + 2.6227x1 という式がそれです。 与えたデータと回帰曲線を同じグラフに表して見ましょう。 これを描くスクリプトはRサンプル4にあります。 dset <- read.csv("d/sample3.csv", header=TRUE) attach (dset) eq1 <- lm(y ~ x1) summary(eq1) plot(x1,y) 縦軸にx1,横軸にyをとって散布図を描く par(new=T) グラフを重ね合わせるため、前のグラフを消さずに次を描く(注 この行はなくても、グラフは重ねて出力されることが分かりました) abline(-6.0930,2.6227) y = -6.0930 + 2.6227x1の直線を引く なお、ablineの引数をeq1の結果から直接取り出すこともできます。スクリプト全体を掲げますが、違っているのは最後の行だけです。 dset <- read.csv("d/sample3.csv", header=TRUE) attach (dset) eq1 <- lm(y ~ x1) names(eq1) summary(eq1) plot(x1,y) par(new=T) abline(eq1$coefficients[1],eq1$coefficients[2]) 残差 †サンプル4を実行した後、R Consoleというウインドウで > z = residuals.lm(eq1) > z と入力してみてください。単に変数名だけ入力すると、その変数の中身が表示されます。 1 2 3 4 5 6 -1.041073 2.060774 -1.673263 1.377661 -2.794107 2.070009 と表示されましたか? データと回帰曲線の(縦方向の)差が残差で、lmの結果につけた名前(eq1 <- lm(y ~ x1)という行があるのがわかりますね?)を使って、上のように取り出すことができます。 t値、P値は何を検定しているか †もしx1の係数が0だったら、 y = -6.0930 という水平線になるはずです。ablineを使って水平線を引くこともできますから、これを同じグラフに表示してみます。ただし水平線があまり端にあると見えにくいので、yの上限15、下限-10を指定します(Rサンプル5)。関係する行はこれ。 plot(x1,y,ylim=c(-10,15)) 出てくるグラフはこんなもののはずです。 もしyとx1の関係が本当はy = -6.0930だとしたら、ここにあるようなデータは出てきそうにありません。x1の係数に関するt値やP値は、そのことを判断するものなのです。 sample5はこんなデータです。x1はsample3と同じですが、yは変化の幅は小さく、しかしはっきり右上がりに作ってあります。 x1 y 1.5 -0.1 2.3 0.3 3.8 0.1 4.2 0.4 5.6 0.6 6.3 0.5 サンプル5で、sample3.csvのかわりにsample5.csvを読み込むように1字書き換えて実行すると、 >Coefficients: > Estimate Std. Error t value Pr(>|t|) >(Intercept) -0.16321 0.16886 -0.967 0.3885 >x1 0.11727 0.03932 2.983 0.0406 * という結果が出てきます。x1の係数は5%有意と判断されています。ついでに水平線もy=-0.16321に置き換えてグラフを出させると、 おっと。グラフの上限下限がさっきのままでした。plot(x1,y,ylim=c(-0.2,0.8))に変更してみましょう(Rサンプル5)。 同じ数字なのに、ふたつのグラフは全然違う印象を与えますね。直観に頼る危険の表れです。 0と0.11727はかなり近いので、x1の係数が0であるという仮説は棄却しにくいのです。しかし残差が小さく、式の当てはまりが非常に良いので、この分散でx1の係数が0だと、こんなにはっきりした右上がりの傾向は出にくいでしょう。そのふたつが影響して、5%有意という少し曖昧な結果が出ています。 |