R言語の回帰分析の実行方法を紹介していきます。回帰分析の実行方法だけでなく、回帰分析により作ったモデルでテストデータの予測も行っていきます。
実行方法では、回帰分析の結果の見方やその抽出方法を詳しく解説し、テストデータの予測などはプロットを用いて視覚的に分かりやすく解説していきます。
今回扱うプログラミングコードも全て、以下からダウンロード可能です。
テンプレートとして自由に扱ってください。
ポチップ単回帰分析
データセットの用意
今回は、次のrockという油田の石に関するデータセットを回帰分析に用います。
1 2 3 4 5 6 7 8 | > head(rock) area peri shape perm 1 4990 2791.90 0.0903296 6.3 2 7002 3892.60 0.1486220 6.3 3 7558 3930.66 0.1833120 6.3 4 7352 3869.32 0.1170630 6.3 5 7943 3948.54 0.1224170 17.1 6 7979 4010.15 0.1670450 17.1 |
rockはnumeric型のarea、peri、shapeとfactor型のpermという変数で構成されています。
このデータセットを関数plotで描画すると、次のように各変数の相関図をプロットをすることができます。
また、この相関の尺度である相関係数は、相関行列を求める関数cov2corと、関数varを用いることで次のように列挙することが可能です。
1 2 3 4 5 6 7 | > corMatrix <- cov2cor(var(dataset)) > corMatrix area peri shape perm area 1.0000000 0.8225064 -0.1821611 -0.3966370 peri 0.8225064 1.0000000 -0.4331255 -0.7387158 shape -0.1821611 -0.4331255 1.0000000 0.5567208 perm -0.3966370 -0.7387158 0.5567208 1.0000000 |
補足として、分散を計算する関数varの引数にデータフレームなどの複数の変数から成るデータを代入すると、共分散行列を計算することが可能です。この共分散行列をcov2corの引数に代入することにより、共分散行列を対応する相関行列に変換することができます。
上の各相関係数（corMatrixの要素）を見ると、areaとperiの相関係数が0.8225064と高く、相関が強いことが分かります。
このように、回帰分析などの統計解析をおこなう前に、データをプロットしたり変数間の相関係数を求めたりすると、どのような解析を行えば良いのかが明確になります。
今回は、areaとperiの相関が強いので、この変数に関する回帰分析を行っていきます。
実行方法
次のように、関数lmの最初の引数にformula、2つ目の引数にデータセットを代入することで、簡単に単回帰分析を実行できます。
1 2 | #単回帰分析 regression <- lm(peri ~ area, data = dataset) |
上の例では、periをareaで回帰しています。
これを統計学のモデルで表すと次のようになります。\begin{align}y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \ \ i = 1, \ldots, n,\end{align}ここに、\(y_i\)はperiを表す確率変数であり、\(y_i\sim N(\beta0 + \beta_1x_1, \sigma^2),\ \ \varepsilon_i \sim N(0, \sigma^2),\ i = 1, \ldots, n\)。
つまり、\(y_i\ (peri)\)と\(\beta_0 + \beta_1x_i\ (area),\ i = 1, \ldots, n\)との差が最小になるような\(\beta_0\)と\(\beta_1\)を計算しています。
最小二乗解である\(\beta_0\)と\(\beta_1\)については、プロットによる解釈で説明しています。
関数lmで得られた変数regressionの中身をみていきます。regressionを参照すると次のモデルと回帰係数が表示されます。
1 2 3 4 5 6 7 8 | > regression Call: lm(formula = peri ~ area, data = dataset) Coefficients: (Intercept) area -471.4358 0.4388 |
さらに詳細な結果が欲しい場合、関数summaryを用います。
1 | summaryReg <- summary(regression) |
上で与えた変数summaryRegを参照してみましょう。次のようにモデルや回帰係数に加え、残差（Std. Error）、t値（t value）、p値（Pr(>|t|)）などの重要な値があることが分かります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | > summaryReg Call: lm(formula = peri ~ area, data = dataset) Residuals: Min 1Q Median 3Q Max -2306.8 -502.3 122.5 564.5 1291.9 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -471.43579 342.77487 -1.375 0.176 area 0.43875 0.04473 9.808 7.51e-13 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 823.1 on 46 degrees of freedom Multiple R-squared: 0.6765, Adjusted R-squared: 0.6695 F-statistic: 96.2 on 1 and 46 DF, p-value: 7.506e-13 |
さらにCoefficientsの項目の変数areaの行に「***」の記号があります。
その下のSignif. codesにあるように、「***」はp値が0.001未満であることが分かります。したがって、変数areaは有意（\beta_1 = 0でない）ことがいえます。
早速、これらの値を取りだしていきます。各検定結果はリストの要素として格納されているため、次のようにsummaryRegの後ろに$coefficientsを付けくわえることで、抽出可能です。
1 2 3 4 5 | >testResults <- summaryReg$coefficients #検定結果 > testResults Estimate Std. Error t value Pr(>|t|) (Intercept) -471.4357929 342.77487413 -1.375351 1.756839e-01 area 0.4387544 0.04473312 9.808268 7.506269e-13 |
上の検定結果testResultsはデータフレームであるので、各列名を与えることで回帰係数やp値の列を取り出すことができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | > testResults[, "Estimate"] #回帰係数 (Intercept) area -471.4357929 0.4387544 > > testResults[, "Std. Error"] #標準誤差 (Intercept) area 342.77487413 0.04473312 > > testResults[, "t value"] #t値 (Intercept) area -1.375351 9.808268 > > testResults[, "Pr(>|t|)"] #p値 (Intercept) area 1.756839e-01 7.506269e-13 |
また、補足として有意な変数を取得したい場合、次のようにp値が有意水準（0.05）未満である行名を得ます。
1 2 3 | > significantVar <- rownames(testResults)[testResults[, "Pr(>|t|)"] < 0.05] > significantVar [1] "area" |
また、目的変数\(Y_i\)に対する説明変数\(x_i\)の情報量の尺度である決定係数は次のようにして得られます。
1 2 3 4 | >Rsquared <- summaryReg$r.squared #決定係数R^2 >adj_Rsquared <- summaryReg$adj.r.squared #修正済み決定係数R^2 > Rsquared [1] 0.6765168 |
決定係数R^2 = 0.6765168でるため、「モデルの当てはまりはまあまあである」ことが言えます。
解析結果の出力
解析結果をcsvファイルに出力してみましょう。
実行方法で得た検定結果をデータフレームにまとめ、write.csvを用いてcsvファイルに出力していきます。
回帰係数、残差、t値、p値、決定係数、修正済み決定係数を解析結果として出力してみます。
次のように関数cbindを用いて、testResultsと決定係数の列Rsquaredと修正済み決定係数の列adjRsquaredを結合させます。
1 2 3 | testResults <- cbind(testResults, Rsquared = c(Rsquared, rep(NA, nrow(testResults) - 1)), adjRsquared = c(adj_Rsquared, rep(NA, nrow(testResults) - 1))) |
testResultsを参照すると、次の解析結果の表が表示されます。
1 2 3 4 | > testResults Estimate Std. Error t value Pr(>|t|) Rsquared adjRsquared (Intercept) -471.4357929 342.77487413 -1.375351 1.756839e-01 0.6765168 0.6694845 area 0.4387544 0.04473312 9.808268 7.506269e-13 NA NA |
最後に、write.csvを用いて、このデータフレームtestResultsを出力します。
1 | write.csv(testResults, "単回帰分析.csv") |
作業ディレクトリに単回帰分析.csvが出力され、次の画像のように単回帰分析の結果が保存されていることが確認できます。
テストデータを用いた予測
回帰分析のモデルの使用例として、テストデータ等の予測があります。
モデルを作ることによって、未知のデータを予測できるようになります。簡単な実行例をみていきましょう。
テストデータとして、次のデータセットvalidationDataを用います。このデータは、さきほどのデータセットrockに正規乱数を加えたものとなります。
1 2 3 4 5 | randomValues <- rnorm(nrow(dataset)) validationData <- data.frame(area = sd(dataset$area) * randomValues + dataset$area, peri = sd(dataset$peri) * randomValues + dataset$peri, shape = sd(dataset$shape) * randomValues + dataset$shap, perm = dataset$perm) |
validationDataを参照すると、次のような乱数が加わったデータであることが確認できます。
1 2 3 4 5 6 7 8 | > head(validationData) area peri shape perm 1 7374.486 4063.870 0.16451266 6.3 2 6096.386 3409.513 0.12044770 6.3 3 4791.727 2455.031 0.09725129 6.3 4 8508.291 4486.127 0.15303604 6.3 5 8503.676 4247.625 0.13986003 17.1 6 3823.757 1793.595 0.03777244 17.1 |
実行方法で得たモデルregressionに、このテストデータの数値をあてはめたい場合、関数predictを用います。
1 | predictData <- predict(regression, validationData) |
最初の引数にlmで得られたオブジェクト、2つ目にあてはめたいデータを代入することで予測することが可能です。
また、次のように予測結果をプロットすることができます。
1 2 3 4 5 | plot(validationData$area, validationData$peri, xlab = "", ylab = "", xlim = c(-2000, 15000), ylim = c(-2000, 5000)) par(new = TRUE) dataOrder <- order(validationData$area) plot(validationData$area[dataOrder], predictData[dataOrder], col = "blue", xlim = c(-2000, 15000), ylim = c(-2000, 5000), xlab = "area(test)", ylab = "peri(test)", axes=FALSE, type = "l", lwd = 1, lty = 1, pch = 1) |
プロット画面に次のような回帰直線が描画されます。青線がテストデータをモデルにあてはめた際の値を通る回帰直線で、黒い点がテストデータの実際の値になります。テストデータであるperiをある程度予測できていることが確認できました。
プロットによる解析結果の解釈
単回帰分析の幾何学的解釈をしていきます。
まず、install.packagesでggplot2をインストールしましょう。このパッケージを用いることで各種統計解析の結果を簡単にプロットすることができます。
1 2 | install.packages("ggplot2") library(ggplot2) |
インストールしライブラリに追加したら、次のように関数ggplotを実行しましょう。
1 2 3 | ggplot(dataset, aes(x = area, y = peri)) + geom_point() + stat_smooth(method = "lm", col = "red") |
変数periとareaのデータのプロット及び回帰直線とその信頼区間の描画をすることができます。上を実行すると次の画像がプロットされます。
dataset$areaの各値\(x_i\)を\(y_i=\beta_0 + \beta_1 x_i = -471.4357929+ 0.4387544\times x_i ,\ i= 1,\ldots, n\)にあてはめた値を通る直線が上の赤い直線です。
また、濃いグレーの領域は\(y_i\)の95%信頼区間を表します。
ポチップ重回帰分析
単回帰分析に続いて、R言語の重回帰分析の実行方法や解析結果の見方を紹介していきます。
単回帰分析のデータセットの用意で用いたデータセットrockを今回も使用します。
重回帰分析により、peri以外の全ての変数でperiを予測することが可能となります。
単回帰分析ではperiをareaで予測しましたが、他の変数を取り入れることでareaだけでは説明できない情報が加わり、予測精度が向上することができます（必ずしも向上しません）。
実行方法
重回帰分析は関数lmを用いることで実行できます。lmの2つ目の引数dataにdatasetを代入し、1つ目の引数のformulaにperi ~.を代入することで、peri以外の変数を説明変数とする重回帰分析を実行できます。
1 | regression <- lm(peri ~., data = dataset) |
これは、次の統計学のモデルを意味します。\begin{align}Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \cdots + \beta_px_{ip} + \varepsilon_i = \boldsymbol{x}_i^T\boldsymbol{\beta},\ \ i = 1, \ldots, n,\end{align}ここに\begin{align}\boldsymbol{x}_i &= \begin{pmatrix}1\\x_{i1}\\ \vdots \\ x_{ip}\end{pmatrix}\ \ \ \ \boldsymbol{\beta} = \begin{pmatrix}\beta_0\\\beta_1\\\vdots\\\beta_p\end{pmatrix}\end{align}であり、\(Y_i\sim N(\boldsymbol{x}_i^T\boldsymbol{\beta}, \sigma^2),\ \ \varepsilon_i \sim N(0, \sigma^2)\)である。
lmの第1引数のformulaを書き換えることで、いろいろなモデルを構築することが可能です。
すなわち、次の実行結果は上と同じです。
1 | regression <- lm(peri ~ area + shape + perm, data = dataset) 上と同じ |
また、+演算子を用いることで次のように任意の説明変数を選択することも可能です。
1 | regression <- lm(peri ~ area + shape, data = dataset) |
これは、areaとshapeの2変数を説明変数とした重回帰分析となります。
また、補足として交互作用ありのモデルは次のように表現できます。
1 | regression <- lm(peri ~ area * shape * perm , data = dataset) 交互作用あり |
*演算子用いることで変数間の交互作用を含めることができます。
関数summaryを用いることで、解析結果得ることができます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | > summaryReg <- summary(regression) > summaryReg Call: lm(formula = peri ~ ., data = dataset) Residuals: Min 1Q Median 3Q Max -1826.58 -211.17 64.82 316.10 1023.32 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.125e+03 3.190e+02 3.525 0.001 ** area 3.368e-01 3.019e-02 11.156 2.06e-14 *** shape -1.150e+03 1.072e+03 -1.072 0.290 perm -1.475e+00 2.191e-01 -6.731 2.84e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 509.3 on 44 degrees of freedom Multiple R-squared: 0.8815, Adjusted R-squared: 0.8734 F-statistic: 109.1 on 3 and 44 DF, p-value: < 2.2e-16 |
上の検定結果の表（Coefficients）の項目を見ると、(Intercept)とareaとpermに「**、***」の記号がついているのが分かります。
その下のSignif. codesの項目を見ると、「0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1」書いており、「**」はp値が0.01未満であり、「***」は0.001未満であることが分かります。
よって、(Intercept)とareaとpermは有意である（\beta_0 = 0、\beta_2 = 0、\beta_3 = 0でない）ことがいえます。
上の解析結果を取り出すには、単回帰分析と同様にリストsummaryRegの各要素にアクセスする必要があります。
次のようにsummaryRegの後ろに$coefficientsを付けることで、回帰係数、残差、t値、p値を取得することができます。
1 2 3 4 5 6 7 | > testResults <- summaryReg$coefficients #検定結果 > testResults Estimate Std. Error t value Pr(>|t|) (Intercept) 1124.5591482 3.190360e+02 3.524866 1.002750e-03 area 0.3368297 3.019366e-02 11.155645 2.060118e-14 shape -1149.8260964 1.072489e+03 -1.072110 2.895160e-01 perm -1.4745435 2.190827e-01 -6.730533 2.841795e-08 |
p値が0.05未満の変数を取得したい場合、次のように"Pr(>|t|)"列が0.05未満であるtestResultsの行名を抽出します。
1 2 3 | > significantVar <- rownames(testResults)[testResults[, "Pr(>|t|)"] < 0.05] > significantVar [1] "(Intercept)" "area" "perm" |
また、目的変数に対する説明変数の情報量の指標である重相関係数は次のようにして得られます。
1 2 3 4 | > Rsquared <- summaryReg$r.squared #重相関係数R^2 > adj_Rsquared <- summaryReg$adj.r.squared #修正済み重相関係数R^2 > Rsquared [1] 0.8815147 |
上の結果から、重相関係数R^2 = 0.8815147であり、モデルの当てはまりが非常に良いことが言えます。
解析結果の出力
解析結果をcsvファイルへ出力する例を紹介します。
単回帰分析と同様に、次のようなCoefficientsの表と重相関係数からなるデータフレームを出力する方法を考えます。
1 2 3 | testResults <- cbind(testResults, Rsquared = c(Rsquared, rep(NA, nrow(testResults) - 1)), adjRsquared = c(adj_Rsquared, rep(NA, nrow(testResults) - 1))) |
データフレームtestResultsには、切片項と各説明変数に対する回帰係数などが含まれているのが分かります。
1 2 3 4 5 6 | > testResults Estimate Std. Error t value Pr(>|t|) Rsquared adjRsquared (Intercept) 1124.5591482 3.190360e+02 3.524866 1.002750e-03 0.8815147 0.8734361 area 0.3368297 3.019366e-02 11.155645 2.060118e-14 NA NA shape -1149.8260964 1.072489e+03 -1.072110 2.895160e-01 NA NA perm -1.4745435 2.190827e-01 -6.730533 2.841795e-08 NA NA |
write.csvを用いて、"重回帰分析.csv"というファイルに出力してみましょう。
1 | write.csv(testResults, "重回帰分析.csv") |
作業ディレクトリを確認すると、"重回帰分析.csv"が出力されているのが分かります。ファイルを開くと解析結果が保存されていることがわかります。
テストデータを用いた予測
モデリングが終わったので、テストデータを先ほどのモデルにあてはめてみましょう。
テストデータを用いた予測で用いたテストデータvalidationDataを再び用います。
validationDataを先ほどのモデルにあてはめたい場合、関数predictを用います。
1 2 3 4 | > predictData <- predict(regression, validationData) > head(predictData) 1 2 3 4 5 6 2789.705 2531.165 3505.780 3820.436 4042.012 3866.934 |
predictの1つ目の引数に関数lmで得られたオブジェクト、2つ目の引数にモデルにあてはめたいデータセット（テストデータ）を代入することで、テストデータの予測が可能です。
上のようにhead(predictData)を参照すると、1行目から6のテストデータをモデルにあてはめた結果が返ってきます。
次にこれらのテストデータとその予測値をプロットして、予測の精度を可視化させましょう。
1 2 3 | plot(validationData$peri, xlab = "", ylab = "") par(new = TRUE) plot(predictData, col = "blue", axes = FALSE, pch = 4, ylab = "peri") |
次を実行することで、テストデータ（黒まる）と予測値（青ばつ）のプロットが可能です。次のような画像がプロットされます。
x軸は標本のインデックス（行番号）であり、だいたい同じ領域に予測値が描画されているのが分かります。
プロットによる解析結果の解釈（説明変数2個）
最後に、重回帰分析の解析結果（回帰平面）をプロットで可視化させます。
単回帰分析の時は、回帰直線でモデルの予測値を描画することができましたが、説明変数が2個以上になると解釈が難しくなります。
今回は、例として説明変数が2個の場合を考えます。次のareaとshapeを説明変数とするモデルを考えます。
1 | regression <- lm(peri ~ area + shape, data = dataset) |
次にプロットの準備として、次のパッケージをインストールします。
1 2 | install.packages("scatterplot3d") library(scatterplot3d) |
パッケージscatterplot3dをインストールすることで、3Dの図形を描画することが可能となります。
パッケージをインストールしたら、グラフのx、z軸をそれぞれarea、shapeに設定し、グラフのy軸を予測値（関数predictの値）に設定することで、回帰平面と実際のデータの座標を描画することができます。
1 2 3 4 | predictGrid <- expand.grid(area = dataset$area, shape = dataset$shape) predictGrid$Volume2 <- predict(regression, new = predictGrid) reg3dspace <- scatterplot3d(pred_grid$area, pred_grid$shape, predictGrid$Volume2, color = "skyblue", pch = 1, ylab = "peri", xlab = "area", zlab = "shape") reg3dspace$points3d(trees$Girth, trees$Height, trees$Volume) |
上を実行することで、次のような画像がプロットされます。
回帰平面実行方法で紹介した統計モデルで、目的変数\(Y_i, \ i = 1, \ldots, n\)はp次元空間上のベクトル\(\boldsymbol{x}_i^T\boldsymbol{\beta}\)であることを仮定しました。
この例では、\(Y_i\)は2次元空間上のベクトルで与えられるため、平面上に分布することが分かります（\(\beta_0 = \beta_1 = \beta_2 = 0\)のとき、\(Y_i=0\)の平面）。
まとめ
R言語で回帰分析の実行方法や解析結果の見方や出力の方法をみていきました。
また、応用としてテストデータを用いた予測や回帰分析のプロットなどを行いました。
次回は、一般化線形モデルの実行方法について解説していきます。
代表的な一般化線形回帰であるロジスティック回帰分析の方法を紹介していきます。
