準備
Sageで回帰分析の説明をするために必要な前準備を以下で行います。
# Rのユーティリティ関数を読み込む
attach(DATA+'RUtil.py')
# プロット関数定義
def showPlot(f):
f.show(aspect_ratio=1, figsize=(5), xmin=-1, xmax=1, ymin=-1, ymax=1)
|
|
回帰分析の注意点
大きな母集団から少ないサンプルを抽出して回帰分析を行う場合、以下のような問題点があります。
- モデルがサンプルに合いすぎる（オーバーフィッティング）問題
- 属性を絞り込む場合、最適なモデルの選択問題
オーバーフィッティングの例
サンプルデータ以下の様な関係を持つとします。 $$ f(x, a) = a_0 + a_1 x $$
ここで、$a_0 = -0.3, a_1 = 0.5$とし、一様部分$U(x|-1,1)$から選んだ$x_n$に、目的値$t_n$に 分散0.3の正規分布の誤差を追加します。 $$ t_n = f(x_n, a) + \mathcal{N}(0, 0.3) $$
# データセット f(x, a) = a0 + a1*x <= a0 = -0.3, a1=0.5
a0 = -0.3;
a1 = 0.5;
# f(x)のプロット
x = var('x')
a_plt = plot(a0 + a1*x, [x, -1, 1], rgbcolor ='red')
|
|
母集団として50個のサンプルデータを作成して、その結果を$f(x,a)$の直線と一緒にプロットしてみます。
データ数が50個の場合には、データにフィットする直線を人の直感で引いた場合とオリジナルの関数$f(x,a)$（赤の直線）はほぼ一致します。
# 観測値は、tn = f(x,a) + N(0, 0.3)
#xv = [];
#tv = [];
#for xn in range(50):
# xn = 2*random() - 1;
# xv += [xn];
# tv += [a0 + a1*xn + gauss(0, 0.3)]
#save(xv, DATA+'xv')
#save(tv, DATA+'tv')
xv = load(DATA+'xv')
tv = load(DATA+'tv')
data = zip(xv, tv);
data_plt = list_plot(zip(xv, tv));
showPlot(a_plt + data_plt)
|
|
50個の母集団から5個のサンプルを抽出すた場合、オリジナルの関数$f(x,a)$からはずれた分布になります。
# 5点を抽出して
data = zip(xv[10:15],tv[10:15])
showPlot(a_plt + list_plot(data))
|
|
<H3>少ないサンプルから求めた最適な直線</H3>
試しに、少ないサンプルから求めた最適な直線を求めています（緑の直線）。
少ないサンプルから求めた最適な直線はオリジナルの関数からずれたり、異なった形になるため、サンプルに合いすぎるのを防ぐ必要があります。
このように回帰分析の結果がサンプルに合いすぎることをオーバーフィッティングと呼びます。
# 5点のサンプルから最適直線を求める
(a0,a1) = var('a0 a1')
model(x) = a0 + a1*x
fit = find_fit(data, model, solution_dict=True); view(fit)
\newcommand{\Bold}[1]{\mathbf{#1}}\left\{a_{0} : -0.130440448397, a_{1} : 0.537229455444\right\}
\newcommand{\Bold}[1]{\mathbf{#1}}\left\{a_{0} : -0.130440448397, a_{1} : 0.537229455444\right\}
|
# 結果の表示
f_fit(x) = model.subs(fit)
fit_plt = plot(f_fit, [x,-1,1], rgbcolor='green')
showPlot(fit_plt + a_plt + list_plot(data))
|
|
Rの線形回帰パッケージ
Rの線形回帰パッケージとして、lmパッケージを使用します。
lmでは、回帰分析の対象の変数を目的変数と呼び、回帰に使用する変数を説明変数と呼んでいます。
lmの最初の引数には、目的変数と説明変数の関係を以下の様な形式で指定します。
目的変数 ~ 説明変数1 + 説明変数2 + ... + 説明変数k
lmを使った直線回帰の例
直線回帰のテストデータとして、女性の身長と体重のサンプルデータwomenを使用して直線回帰を実施してみます。
回帰の結果は、fit変数に保存され、summary関数で回帰の結果を確認できます。
回帰の結果求まった係数（Estimate列）がどの程度重要なのかは、行の最後の重要度（***）の部分で判断できます。★が多いほど重要であることを示しています。
このようにlmパッケージを使った線形回帰分析では、係数の算出だけでなく、その重要度までも含めた情報を出力します。
# Rの線形回帰パッケージを使った分析
r('fit <- lm(weight ~ height, data=women)')
r('summary(fit)')
Call:
lm(formula = weight ~ height, data = women)
Residuals:
Min 1Q Median 3Q Max
-1.7333 -1.1333 -0.3833 0.7417 3.1167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
height 3.45000 0.09114 37.85 1.09e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14
Call:
lm(formula = weight ~ height, data = women)
Residuals:
Min 1Q Median 3Q Max
-1.7333 -1.1333 -0.3833 0.7417 3.1167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
height 3.45000 0.09114 37.85 1.09e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14
|
graph = preGraph("fig8_1.pdf")
r('plot(women$height, women$weight,xlab="Height (in inches)",ylab="Weight (in pounds)")')
r('abline(fit)')
postGraph(graph, fac=0.5)
|
|
# 多重回帰分析の例
# r('install.packages("car")')
r('library(car)')
[1] "car" "nnet" "MASS" "stats" "graphics" "grDevices" "utils" "datasets" [9] "methods" "base" [1] "car" "nnet" "MASS" "stats" "graphics" "grDevices" "utils" "datasets" [9] "methods" "base" |
分析用のデータ
多重回帰分析のテストデータとして、Rに付属するアメリカ合衆国のstate.x77の中から、犯罪率（Murder）、人口（Population）、識字率（Illiteracy）、平均収入（Income）、真冬日の日数（Frost）を取り出し、犯罪率との関係を回帰分析で求めます。
# データを用意する
r('states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])')
Murder Population Illiteracy Income Frost Alabama 15.1 3615 2.1 3624 20 Alaska 11.3 365 1.5 6315 152 Arizona 7.8 2212 1.8 4530 15 Arkansas 10.1 2110 1.9 3378 65 California 10.3 21198 1.1 5114 20 Colorado 6.8 2541 0.7 4884 166 Connecticut 3.1 3100 1.1 5348 139 Delaware 6.2 579 0.9 4809 103 Florida 10.7 8277 1.3 4815 11 Georgia 13.9 4931 2.0 4091 60 Hawaii 6.2 868 1.9 4963 0 Idaho 5.3 813 0.6 4119 126 Illinois 10.3 11197 0.9 5107 127 Indiana 7.1 5313 0.7 4458 122 Iowa 2.3 2861 0.5 4628 140 Kansas 4.5 2280 0.6 4669 114 Kentucky 10.6 3387 1.6 3712 95 Louisiana 13.2 3806 2.8 3545 12 Maine 2.7 1058 0.7 3694 161 Maryland 8.5 4122 0.9 5299 101 Massachusetts 3.3 5814 1.1 4755 103 Michigan 11.1 9111 0.9 4751 125 Minnesota 2.3 3921 0.6 4675 160 Mississippi 12.5 2341 2.4 3098 50 Missouri 9.3 4767 0.8 4254 108 Montana 5.0 746 0.6 4347 155 Nebraska 2.9 1544 0.6 4508 139 Nevada 11.5 590 0.5 5149 188 New Hampshire 3.3 812 0.7 4281 174 New Jersey 5.2 7333 1.1 5237 115 New Mexico 9.7 1144 2.2 3601 120 New York 10.9 18076 1.4 4903 82 North Carolina 11.1 5441 1.8 3875 80 North Dakota 1.4 637 0.8 5087 186 Ohio 7.4 10735 0.8 4561 124 Oklahoma 6.4 2715 1.1 3983 82 Oregon 4.2 2284 0.6 4660 44 Pennsylvania 6.1 11860 1.0 4449 126 Rhode Island 2.4 931 1.3 4558 127 South Carolina 11.6 2816 2.3 3635 65 South Dakota 1.7 681 0.5 4167 172 Tennessee 11.0 4173 1.7 3821 70 Texas 12.2 12237 2.2 4188 35 Utah 4.5 1203 0.6 4022 137 Vermont 5.5 472 0.6 3907 168 Virginia 9.5 4981 1.4 4701 85 Washington 4.3 3559 0.6 4864 32 West Virginia 6.7 1799 1.4 3617 100 Wisconsin 3.0 4589 0.7 4468 149 Wyoming 6.9 376 0.6 4566 173 Murder Population Illiteracy Income Frost Alabama 15.1 3615 2.1 3624 20 Alaska 11.3 365 1.5 6315 152 Arizona 7.8 2212 1.8 4530 15 Arkansas 10.1 2110 1.9 3378 65 California 10.3 21198 1.1 5114 20 Colorado 6.8 2541 0.7 4884 166 Connecticut 3.1 3100 1.1 5348 139 Delaware 6.2 579 0.9 4809 103 Florida 10.7 8277 1.3 4815 11 Georgia 13.9 4931 2.0 4091 60 Hawaii 6.2 868 1.9 4963 0 Idaho 5.3 813 0.6 4119 126 Illinois 10.3 11197 0.9 5107 127 Indiana 7.1 5313 0.7 4458 122 Iowa 2.3 2861 0.5 4628 140 Kansas 4.5 2280 0.6 4669 114 Kentucky 10.6 3387 1.6 3712 95 Louisiana 13.2 3806 2.8 3545 12 Maine 2.7 1058 0.7 3694 161 Maryland 8.5 4122 0.9 5299 101 Massachusetts 3.3 5814 1.1 4755 103 Michigan 11.1 9111 0.9 4751 125 Minnesota 2.3 3921 0.6 4675 160 Mississippi 12.5 2341 2.4 3098 50 Missouri 9.3 4767 0.8 4254 108 Montana 5.0 746 0.6 4347 155 Nebraska 2.9 1544 0.6 4508 139 Nevada 11.5 590 0.5 5149 188 New Hampshire 3.3 812 0.7 4281 174 New Jersey 5.2 7333 1.1 5237 115 New Mexico 9.7 1144 2.2 3601 120 New York 10.9 18076 1.4 4903 82 North Carolina 11.1 5441 1.8 3875 80 North Dakota 1.4 637 0.8 5087 186 Ohio 7.4 10735 0.8 4561 124 Oklahoma 6.4 2715 1.1 3983 82 Oregon 4.2 2284 0.6 4660 44 Pennsylvania 6.1 11860 1.0 4449 126 Rhode Island 2.4 931 1.3 4558 127 South Carolina 11.6 2816 2.3 3635 65 South Dakota 1.7 681 0.5 4167 172 Tennessee 11.0 4173 1.7 3821 70 Texas 12.2 12237 2.2 4188 35 Utah 4.5 1203 0.6 4022 137 Vermont 5.5 472 0.6 3907 168 Virginia 9.5 4981 1.4 4701 85 Washington 4.3 3559 0.6 4864 32 West Virginia 6.7 1799 1.4 3617 100 Wisconsin 3.0 4589 0.7 4468 149 Wyoming 6.9 376 0.6 4566 173 |
分布図で各変数の関連を調べる
carライブラリのscatterplotMatrix関数を使って各変数の関連を調べます。 MurderとIlliteracyに正の相関が見られ、後は判断が難しいように思われます。
# 分布図で関係を見る
graph = preGraph("fig8_4.pdf")
r('scatterplotMatrix(states, spread=FALSE, lty.smooth=2)')
postGraph(graph, fac=0.7)
|
|
多重回帰の結果
lm関数を使って上記のデータの多重回帰分析を行います。
処理は、ほとんど線形回帰と同じです。解析結果のsummaryも分布図から見て取られる傾向と合致し、犯罪率（Murder）には、識字率（Illiteracy）と人口（Population）が関連があることが分かります。
# lmを使った多重回帰
r('fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=states)')
r('summary(fit)')
Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
data = states)
Residuals:
Min 1Q Median 3Q Max
-4.7960 -1.6495 -0.0811 1.4815 7.6210
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.235e+00 3.866e+00 0.319 0.7510
Population 2.237e-04 9.052e-05 2.471 0.0173 *
Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
Income 6.442e-05 6.837e-04 0.094 0.9253
Frost 5.813e-04 1.005e-02 0.058 0.9541
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
data = states)
Residuals:
Min 1Q Median 3Q Max
-4.7960 -1.6495 -0.0811 1.4815 7.6210
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.235e+00 3.866e+00 0.319 0.7510
Population 2.237e-04 9.052e-05 2.471 0.0173 *
Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
Income 6.442e-05 6.837e-04 0.094 0.9253
Frost 5.813e-04 1.005e-02 0.058 0.9541
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
|
解析結果のグラフ
解析結果fitをplotすると、予測値と実際の値との残差の関係（左上）、残差のQ-Qプロット（右上）、残差の平方根プロット（左下）、残差と影響力プロット（右下）が表示される。
- 線型モデルがデータを表現できていれば、残差（Residuals）と予測値（Fitted values）はランダムに分布する
- 線型モデルの仮定である残差が正規分布かどうかは、Q-Qプロットが45度の直線に乗っているかどうかで判断する
- 分散を定数だと仮定すると、残差の平方根プロットは、水平線にランダムにプロットされる
- データの中で特に逸脱している（影響を及ぼしている）データをクック距離（Cook's distanceの破線）から推定する。クック距離が0.5以上だと影響力が大きい（例ではAlaska, Nevadaのデータが大きい）
# 回帰結果のプロット例
graph = preGraph("fig8_7.pdf")
r('par(mfrow=c(2,2))')
r('plot(fit)')
postGraph(graph, fac=0.7)
|
|
分析結果の評価
回帰分析の結果から各説明変数が1%変化した場合、目的変数の変化する範囲を95%の確率で求めた結果がconfint関数の結果です。
例えば識字率（Illiteracy）を1%変化させると犯罪率（Murder）は、2.38から5.90の範囲で変化すると予測されます。
# 分析結果の評価
r('confint(fit)')
2.5 % 97.5 % (Intercept) -6.552191e+00 9.0213182149 Population 4.136397e-05 0.0004059867 Illiteracy 2.381799e+00 5.9038743192 Income -1.312611e-03 0.0014414600 Frost -1.966781e-02 0.0208304170 2.5 % 97.5 % (Intercept) -6.552191e+00 9.0213182149 Population 4.136397e-05 0.0004059867 Illiteracy 2.381799e+00 5.9038743192 Income -1.312611e-03 0.0014414600 Frost -1.966781e-02 0.0208304170 |
最適なモデルの選別
線型モデルの最適なモデルの選択を行うために、赤池情報量基準（AIC）を使って線型モデルの説明変数の組み合わせから、最適なモデルを見つける関数がMASSパッケージのstepAIC関数です。
stepAICの結果、PopulationとIlliteracyのモデルが最適なモデルであることが求まりました。
# AICをベースに最適なモデルを求める
r('library(MASS)')
r('fit<-lm(Murder ~ Population + Illiteracy + Income + Frost,data=states)')
r('stepAIC(fit, direction="backward")')
Call: lm(formula = Murder ~ Population + Illiteracy, data = states) Coefficients: (Intercept) Population Illiteracy 1.6515497 0.0002242 4.0807366 Call: lm(formula = Murder ~ Population + Illiteracy, data = states) Coefficients: (Intercept) Population Illiteracy 1.6515497 0.0002242 4.0807366 |
K-分割交差検証
オーバーフィッティングを防ぐためにK-分割交差検証を使用します。
K-分割交差検証は、サンプルデータをk個に分割し、K-1個を訓練データに、残る1個を検証用データとして回帰分析を繰り返し、K回の結果を平均して解を求めます。
#r('install.packages("bootstrap")')
r('library(bootstrap)')
# shrinkage関数をsourceコマンドで読み込む
filename = DATA+"shrinkage.txt"
r('source("%s")' % filename)
$value
function (fit, k = 10)
{
theta.fit <- function(x, y) {
lsfit(x, y)
}
theta.predict <- function(fit, x) {
cbind(1, x) %*% fit$coef
}
x <- fit$model[, 2:ncol(fit$model)]
y <- fit$model[, 1]
results <- crossval(x, y, theta.fit, theta.predict, ngroup = k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
ret <- paste("Original R-square =", r2, "\n")
ret <- paste(ret, k, "Fold Cross-Validated R-square =", r2cv,
"\n")
ret <- paste(ret, "Change =", r2 - r2cv, "\n")
return(ret)
}
$visible
[1] FALSE
$value
function (fit, k = 10)
{
theta.fit <- function(x, y) {
lsfit(x, y)
}
theta.predict <- function(fit, x) {
cbind(1, x) %*% fit$coef
}
x <- fit$model[, 2:ncol(fit$model)]
y <- fit$model[, 1]
results <- crossval(x, y, theta.fit, theta.predict, ngroup = k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
ret <- paste("Original R-square =", r2, "\n")
ret <- paste(ret, k, "Fold Cross-Validated R-square =", r2cv,
"\n")
ret <- paste(ret, "Change =", r2 - r2cv, "\n")
return(ret)
}
$visible
[1] FALSE
|
例では、shrinkage関数内でK-分割交差検証用の関数crossvalを使って、R-squareが0.06改善されたことを示しています。
# AICで求めたモデルに対し回帰分析を実施
r('fit <- lm(Murder ~ Population + Illiteracy, data=states)')
# shrinkageの実行
r('shrinkage(fit)')
[1] "Original R-square = 0.566832677066343 \n 10 Fold Cross-Validated R-square = 0.541458786180219 \n Change = 0.0253738908861243 \n" [1] "Original R-square = 0.566832677066343 \n 10 Fold Cross-Validated R-square = 0.541458786180219 \n Change = 0.0253738908861243 \n" |
|
|