d01 <- read.csv("practice.csv")『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第10章 重回帰分析（Ⅱ）
0 全体の章構成
1 重回帰モデルの複雑化
1-1 二乗項(squared term)
\[
y=\alpha_{0}+\alpha_{1}(x_{1}-\gamma)^2=(\alpha_{0}+\alpha_{1}\gamma^2)-2\alpha_{1}\gamma\;x_{1}+\alpha_{1}x_{1}^2
\] \((\alpha_{0}+\alpha_{1}\gamma^2)=\beta_{0}\)，\(-2\alpha_{1}\gamma=\beta_{1}\)，\(x_{1}^2=x_{2}\)，\(\alpha_{1}=\beta_{2}\)と書き換えれば， \[y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}\] と云う普通の重回帰の式になる1。
統計分析ソフトはこの後者の様な重回帰モデルとして係数を計算するだけである2。
計算結果から放物線の頂点を逆算すると以下のようになる。 \[(\gamma,\alpha_{0})=(-\frac{\beta_{1}}{2\alpha_{1}},\beta_{0}-\alpha_{1}\gamma^2)=(-\frac{\beta_{1}}{2\beta_{2}},\beta_{0}-\alpha_{1}\frac{\beta_{1}^2}{4\alpha_{1}^2})=(-\frac{\beta_{1}}{2\beta_{2}},\beta_{0}-\frac{\beta_{1}^2}{4\beta_{2}})\] 本文の例を，サポートウェブの表記や流儀に変えて実行してみよう3。
sum100101 <- summary(lm(income ~ age + edu + I(edu^2), data = d01))
sum100101
Call:
lm(formula = income ~ age + edu + I(edu^2), data = d01)
Residuals:
Min 1Q Median 3Q Max
-546.11 -227.55 -59.18 171.87 1958.94
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1466.440 631.237 2.323 0.020773 *
age 5.053 2.243 2.252 0.024941 *
edu -258.732 92.445 -2.799 0.005428 **
I(edu^2) 11.411 3.392 3.365 0.000856 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 320.4 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.1173, Adjusted R-squared: 0.1094
F-statistic: 14.8 on 3 and 334 DF, p-value: 4.537e-09教育年数eduの1次の項も二次の項も1%水準で有意なので，教育年数は個人年収に対して曲線的な関係を有していると考えられる。
予測値のグラフ
with(d01, table(edu)) # x軸の取り得る値を確認edu
9 11 12 14 16 18
15 3 110 123 113 12 .edu <- c(9:18) # x軸の描画範囲を指定
.y <- coef100101[1] +
coef100101[2] * mean(d01$age, na.rm = T) + # 全体の平均年齢を代入
coef100101[3] * .edu +
coef100101[4] * (.edu^2)
RorB <- c("red", "grey", "red", "red", "grey", "red", "grey", "red", "grey", "red")
Mark <- c(15, 0, 15, 15, 0, 15, 0, 15, 0, 15)
plot(.edu, .y,
xlim = c(range(d01$edu, na.rm = T)),
pch = Mark,
col = RorB,
bty = "n",
xlab = "教育年数",
ylab = "収入",
ylim = c(200, 800),
las = 1,
main = "教育年数と年収の予測値の関係（全体の平均年齢で計算）")
単純に直線的に上昇しないのが不思議に思われるだろうが，データの中に専修学校高等課程（教育年数11年として計算）の回答者がわずかに含まれ，この人たちの平均年収が中卒や高卒以上の人達と比べてかなり低いためにこうした結果になっている。
ケース数は上で確認したが、教育年数カテゴリ別平均年収を確認してみよう。
(.observed <- with(d01, by(income, edu, mean, na.rm = T)))edu: 9
[1] 273.3333
------------------------------------------------------------
edu: 11
[1] 208.3333
------------------------------------------------------------
edu: 12
[1] 287.3786
------------------------------------------------------------
edu: 14
[1] 236.4865
------------------------------------------------------------
edu: 16
[1] 537.5
------------------------------------------------------------
edu: 18
[1] 530因みに，それぞれの教育年数グループで平均年齢が異なるので、全体での平均年齢ではなくグループごとの平均年齢を用いて予測値を求める方がより適切だと考えられる。
グループ別平均年齢
# 教育年数のそれぞれのカテゴリごとの平均年齢を用いて年収の予測値を計算する
v_edu <- c(9, 11 ,12, 14, 16, 18)
m_age <- with(d01, by(age, edu, mean, na.rm = T))
m_age # 教育年数グループごとの平均年齢edu: 9
[1] 50.53333
------------------------------------------------------------
edu: 11
[1] 41
------------------------------------------------------------
edu: 12
[1] 47.24545
------------------------------------------------------------
edu: 14
[1] 45.94309
------------------------------------------------------------
edu: 16
[1] 45.0354
------------------------------------------------------------
edu: 18
[1] 36.33333グループ別平均年齢を用いた予測値
.predict <- coef100101[1] +
coef100101[2] * m_age +
coef100101[3] * v_edu +
coef100101[4] * v_edu^2
.predictedu: 9
[1] 317.4867
------------------------------------------------------------
edu: 11
[1] 208.3056
------------------------------------------------------------
edu: 12
[1] 243.5885
------------------------------------------------------------
edu: 14
[1] 312.9305
------------------------------------------------------------
edu: 16
[1] 475.557
------------------------------------------------------------
edu: 18
[1] 690.092この予測値ととグループごとの実測値平均を重ねると以下の通りである。
3つ重ね描きしたグラフ
plot(v_edu, .predict,
xlim = c(range(d01$edu, na.rm = T)),
ylim = c(200, 800),
pch = 15,
col = "#ff000080",
bty = "n",
xlab = "教育年数",
ylab = "収入",
las = 1,
main = paste("年収の予測値（赤■，群別平均年齢で計算）と実測値（青●）\n（赤の放物線は全体平均年齢で作図）"))
par(new = T) # グラフを重ねる為の指定
plot(.x0 <- c(90:180)/10,
coef100101[1] + coef100101[2] * mean(d01$age, na.rm = T) + coef100101[3] * .x0 + coef100101[4] * .x0^2,
xlim = c(range(d01$edu, na.rm = T)),
ylim = c(200, 800),
bty = "n",
type = "l",
axes = F,
col = "#ff000060",
xlab = "",
ylab = "",
main = "")
par(new = T)
plot(v_edu, .observed,
xlim = c(range(d01$edu, na.rm = T)),
ylim = c(200, 800),
bty = "n",
axes = F,
pch = 16,
col = "#0000ff80",
xlab = "",
ylab = "",
main = "")
1-2 交互作用項(interaction term)
\[\begin{alignat*}{3} y_{i} &=\beta_{0}+\beta_{1}\color{red}d_{i}\color{black}&&+\beta_{2}x_{2i}+\beta_{3}\color{red}d_{i}\color{black}x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i},\;\;\;\;(\color{red}d_{i}\color{black}=0\;or\;1)\\ \color{red}d_{i}\color{black}=0\Rightarrow y_{i} &=\beta_{0}&&+\beta_{2}x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i}\\ \color{red}d_{i}\color{black}=1\Rightarrow y_{i} &=\beta_{0}+\beta_{1}&&+\beta_{2}x_{2i}+\beta_{3}x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i}\\ &=\beta_{0}+\beta_{1}&&+(\beta_{2}+\beta_{3})x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i}\\ \end{alignat*}\]
本文では，男性ダミーの主効果を投入して個人年収を説明する重回帰分析で，教育年数は二乗項のみを入れている。その理由も含めて以下に結果を示す。
教育年数の一次の項と二次の項を同時投入する4といずれも有意にならず，分散説明率は.3467（自由度調整済みで.3389）である。
教育年数の一次の項のみだと.3465(.3406)，二次の項のみだと.3467(.3408)。
よって，僅かではあるが二乗項のみの投入が最善である。
d01$male <- c(1, 0)[d01$sex]summary(lm(income ~ male + age + edu + I(edu^2), data = d01))
Call:
lm(formula = income ~ male + age + edu + I(edu^2), data = d01)
Residuals:
Min 1Q Median 3Q Max
-718.02 -169.48 -46.59 112.37 1787.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -421.531 571.202 -0.738 0.46105
male 349.208 32.293 10.814 < 2e-16 ***
age 5.393 1.933 2.790 0.00557 **
edu 11.373 83.474 0.136 0.89171
I(edu^2) 1.134 3.073 0.369 0.71246
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 276.1 on 333 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3467, Adjusted R-squared: 0.3389
F-statistic: 44.19 on 4 and 333 DF, p-value: < 2.2e-16summary(lm(income ~ male + age + edu, data = d01))
Call:
lm(formula = income ~ male + age + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -624.486 153.254 -4.075 5.76e-05 ***
male 352.892 30.670 11.506 < 2e-16 ***
age 5.366 1.929 2.782 0.00571 **
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16summary(reg100102a <- lm(income ~ male + age + I(edu^2), data = d01))
Call:
lm(formula = income ~ male + age + I(edu^2), data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.95 -169.43 -47.52 108.84 1787.45
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -345.3163 115.2701 -2.996 0.00294 **
male 347.8911 30.7684 11.307 < 2e-16 ***
age 5.3980 1.9297 2.797 0.00545 **
I(edu^2) 1.5503 0.2847 5.446 1e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3467, Adjusted R-squared: 0.3408
F-statistic: 59.08 on 3 and 334 DF, p-value: < 2.2e-16この結果から幾つかのパタンについて予測値を計算してみよう。
これについて上と同じ要領で予測値を求めてみる。
ただし，定数項，男性ダミー主効果，年齢主効果は10%水準でも有意ではない（＝0でないとは言えない）ので，計算には含めない。
変数のセンタリング
ここでは，結果を読み取り易くする為に，0を取り得ない値を，「平均からの偏差」に変換して分析してみよう。
現在の結果では定数項も大幅にマイナスとなっており，意味のある値を表現していない。
また，交互作用項の解釈が平易になる様に，男性ダミーではなく女性ダミーとして投入しよう。
\[\begin{alignat*}{3} &\color{red}{no-centering}&&&&\\ y_{i} &=\beta_{0}+\beta_{1}age_{i}&&+\beta_{2}female_{i}+\beta_{3}edu_{i}&&+\beta_{4}female_{i}\times edu_{i}&&+\epsilon_{i}\\ &\color{green}{centering}&&&&\\ y_{i} &=\beta_{0}'+\beta_{1}(age_{i}-\overline{age})&&+\beta_{2}female_{i}+\beta_{3}(edu_{i}-\overline{edu})&&+\beta_{4}female_{i}\times (edu_{i}-\overline{edu})&&+\epsilon_{i}\\ \end{alignat*}\]
年齢の平均
mean(d01$age, na.rm = T)[1] 45.80208d01$d_age <- d01$age - mean(d01$age, na.rm = T)教育年数の平均
mean(d01$edu, na.rm = T)[1] 13.92021d01$d_edu <- d01$edu - mean(d01$edu, na.rm = T)女性ダミー
d01$female <- c(0, 1)[d01$sex]定数項の値は，年齢，教育年数共に平均値に等しい（つまり各偏差が0の）男性（女性ダミーが0）の世帯年収平均を意味する。
性別の主効果は消えてしまった。
本人の教育年数によって世帯年収平均は上昇するが，その上昇幅は女性の方が2倍余りもある（主効果と交互作用効果の比較）。
標準関数による交互作用プロット
標準関数でも，交互作用（従属変数と2つの独立変数の計3変数の関連）を可視化する事が出来る。
以下では，NAを含む場合も問題なく使える様に若干工夫している5。
また，関数の中で数値型変数sexを文字型に変換している6。
標準関数のinteraction.plot
with(d01,
interaction.plot(edu, c("男性", "女性")[sex], fincome,
fun = function(x) mean(x, na.rm = T),
ylim = c(0, 2000),
bty = "l",
lty = c(2, 1),
main = "性・教育年数別の世帯年収平均",
xlab = "本人教育年数（カテゴリカル；間隔が要調整）",
ylab = "平均世帯年収（万円）",
trace.label = c("性別"),
fixed = T,
family = "serif",
las = 1,
col = c("#008b8b", "#deb887"))
)
1-3 多重共線性(multi-collinearity)
\[y=\beta_{0}+\beta_{1}x_{1}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}+\epsilon\] この時， \[x_{j}=\alpha_{0}+\alpha_{1}x_{1}+...+\alpha_{j-1}x_{j-1}+\alpha_{j+1}x_{j+1}+...+\alpha_{p}x_{p}+\delta\] の分散説明率を\(R_{x_{j}}^2\)とする。この時，トレランスとVIF(Variance Inflation Factor)は以下の様に定義される。 \[\begin{align*} tolerance_{x_{j}}&=1-R_{x_{j}}^2\\ VIF_{x_{j}}&=\frac{1}{tolerance_{x_{j}}}=\frac{1}{1-R_{x_{j}}^2}\\ \end{align*}\]
トレランスは小さい程，VIFは大きい程問題とされる。目安は明確に定められないが，経験的には，\(VIF>2\)となったら少し気に留める方が良いかも知れない。
第9章で説明しているが，（標本）偏回帰係数\(\hat{\beta}_{x_{j}}=b_{j}\)の標準誤差の推定値は，下記の式で表される。 \[\hat{\sigma}_{b_{j}}=\frac{\hat{\sigma}_{\epsilon}}{\sqrt{SS_{x_j}}\sqrt{1-R_{x_{j}}^2}}=\frac{\sqrt{MS_{residual}}}{\sqrt{SS_{x_j}}\sqrt{1-R_{x_{j}}^2}}=\frac{\sqrt{MS_{residual}}}{\sqrt{SS_{x_j}}}\sqrt{VIF_{x_{j}}}\] \(H_{0}:\;\beta_{x_{j}}=0\)の検定統計量tは， \[t=\frac{b_{j}}{\hat{\sigma}_{b_{j}}}\sim\;t_{(n-p-1)}\] であるが，分母の標準誤差推定値が大きければ絶対値は小さくなりがちであり，偏回帰係数の推定値自体が大きくても統計的に有意にならない事が有り得る。
「点推定値は大きいのに，統計的に有意にならない！」と云うのがこの問題なのである。
点推定値が大きくても，信頼区間が非常に広くなって0をまたいでしまうのである。
本文と同様の重回帰分析を行って，toleranceやvifを計算してみよう。
library(car)summary(reg100103a <- lm(income ~ age + male + edu, data = d01))
Call:
lm(formula = income ~ age + male + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -624.486 153.254 -4.075 5.76e-05 ***
age 5.366 1.929 2.782 0.00571 **
male 352.892 30.670 11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16car::vif(reg100103a) age male edu
1.051380 1.012674 1.060861 1/vif(reg100103a) # これがトレランスの値 age male edu
0.9511306 0.9874847 0.9426309 モデルreg100103bについては，VIFは大きいものの，標準誤差が著しく大きくなってt検定が有意にならないという問題は生じていない7。
モデルreg100103cでは特にmaleの主効果で標準誤差が著しく大きいので，これを独立変数から除外しても良いかも知れない。
定義に従ったVIFの計算
本文で触れた，carパッケイジを使わずにvifを計算する方法を紹介しておく。
本項冒頭のモデルreg100103aについて自力でVIFを求めるには下記の様にする。
ただしこの方法は，特に独立変数の数が多い場合には手間がかかる8。
reg100103a$calllm(formula = income ~ age + male + edu, data = d01)Full <- with(d01, complete.cases(income, age, male, edu))1-4 変数のコントロールとは
本文で例示している予測値を求める。既に何度か紹介した方法である。
sum(coef100103a[, 1] * c(1, 40, 0, 14)) # 40歳、女性、短大卒[1] 178.5939sum(coef100103a[, 1] * c(1, 40, 0, 16)) # 40歳、女性、四大卒[1] 262.6563sum(coef100103a[, 1] * c(1, 40, 1, 16)) # 40歳、男性、四大卒[1] 615.5484 各カテゴリの実際の標本平均は次の様にして求める事が出来る。
ここでは集計対象が少なくなりすぎるので，40歳ちょうどではなく，35歳から44歳を対象として計算してみた。
# cases1，2，3は，各ケースが指定の条件に該当するか否かをTRUE(1)かFALSE(0)で表現するヴェクトル
cases1 <- d01$sex == 2 & d01$age >= 35 & d01$age <= 44 & d01$edu == 14
cases2 <- d01$sex == 2 & d01$age >= 35 & d01$age <= 44 & d01$edu == 16
cases3 <- d01$sex == 1 & d01$age >= 35 & d01$age <= 44 & d01$edu == 16summary(d01$income[cases1]) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 25.0 125.0 149.2 200.0 800.0 2 summary(d01$income[cases2]) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 25.0 200.0 432.7 600.0 1750.0 3 summary(d01$income[cases3]) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 500.0 600.0 627.6 862.5 1125.0 5 発展1 一般線型モデル(GLM, LM)
かつては一般線型モデル\((General\;Linear\; Model)\)の事を，GLMと呼ぶ事が一般的であったが，現在では一般化線型モデル\((General\color{red}{ized}\;Linear\;Model)\)の事をGLMと略す事が多く，一般線型モデルの方は単にLMと呼ぶ事が多くなって来ている。
性別変数を3通りの方法で重回帰分析に投入したのが以下である。それぞれの結果をよく見比べて，何処が違うかを確認した上で，1番目と3番目についてはそれぞれ予測式を具体的に考えて，本質的に同じである事を確認しよう。
summary(lm(income ~ age + sex + edu, data = d01))
Call:
lm(formula = income ~ age + sex + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 81.298 163.342 0.498 0.61901
age 5.366 1.929 2.782 0.00571 **
sex -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16summary(lm(income ~ age + factor(sex) + edu, data = d01))
Call:
lm(formula = income ~ age + factor(sex) + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
factor(sex)2 -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16d01$female <- c(0, 1)[d01$sex]
summary(lm(income ~ age + female + edu, data = d01))
Call:
lm(formula = income ~ age + female + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
female -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16本文で言及している残りのモデルについても結果を掲示しておくので読み解いてみよう。（表示結果を簡略化する為にsummary( )ではなくlm( )の結果だけを示す。）
lm(income ~ age * sex + edu, data = d01)
Call:
lm(formula = income ~ age * sex + edu, data = d01)
Coefficients:
(Intercept) age sex edu age:sex
-497.3846 17.5482 0.1679 43.6989 -7.7188 lm(income ~ age * factor(sex) + edu, data = d01)
Call:
lm(formula = income ~ age * factor(sex) + edu, data = d01)
Coefficients:
(Intercept) age factor(sex)2 edu
-497.2167 9.8293 0.1679 43.6989
age:factor(sex)2
-7.7188 lm(income ~ age + edu + age:sex , data = d01)
Call:
lm(formula = income ~ age + edu + age:sex, data = d01)
Coefficients:
(Intercept) age edu age:sex
-497.108 17.543 43.698 -7.715 lm(income ~ age + edu + age:factor(sex) , data = d01)
Call:
lm(formula = income ~ age + edu + age:factor(sex), data = d01)
Coefficients:
(Intercept) age edu age:factor(sex)2
-497.108 9.827 43.698 -7.715 発展2 ’良い’モデルの検討
発展2-1 変数選択
まずは本文と同様の変数減少法の結果を紹介した後，それ以外の変数選択法についても簡単に紹介する。
出力そのままではstepごとの区切りが分かりにくいかも知れないので，各ステップをピンクの枠線で囲んである。が何もしない場合であり，どれかの独立変数を除去したらAICが低下するかどうかを確かめながら，それ以上除去する独立変数がなくなった時点で終了する。最後に，採用されたモデルの偏回帰係数が出力される。
step(lm(fincome ~ age * edu * factor(sex), data = d01))Start: AIC=3585.14
fincome ~ age * edu * factor(sex)
Df Sum of Sq RSS AIC
- age:edu:factor(sex) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome ~ age + edu + factor(sex) + age:edu + age:factor(sex) +
edu:factor(sex)
Df Sum of Sq RSS AIC
- age:factor(sex) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age:edu 1 393611 55577691 3584.0
- edu:factor(sex) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome ~ age + edu + factor(sex) + age:edu + edu:factor(sex)
Df Sum of Sq RSS AIC
- age:edu 1 371070 55578768 3582.0
<none> 55207698 3582.1
- edu:factor(sex) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome ~ age + edu + factor(sex) + edu:factor(sex)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
- edu:factor(sex) 1 1003253 56582021 3585.3
- age 1 1606902 57185670 3588.4
Call:
lm(formula = fincome ~ age + edu + factor(sex) + edu:factor(sex),
data = d01)
Coefficients:
(Intercept) age edu factor(sex)2
-481.771 9.573 53.000 -790.292
edu:factor(sex)2
61.381 これをst1というオブジェクトに格納してsummary( )を見ると以下の通りである。また，各ステップの要点だけを表形式に表現するには，オブジェクトの中のanova情報を表示させる。オブジェクトst1に格納する際に上記と同じ変数選択の過程が出力されるが以下では省略する。
st1 <- step(lm(fincome ~ age * edu * factor(sex), data = d01))Start: AIC=3585.14
fincome ~ age * edu * factor(sex)
Df Sum of Sq RSS AIC
- age:edu:factor(sex) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome ~ age + edu + factor(sex) + age:edu + age:factor(sex) +
edu:factor(sex)
Df Sum of Sq RSS AIC
- age:factor(sex) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age:edu 1 393611 55577691 3584.0
- edu:factor(sex) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome ~ age + edu + factor(sex) + age:edu + edu:factor(sex)
Df Sum of Sq RSS AIC
- age:edu 1 371070 55578768 3582.0
<none> 55207698 3582.1
- edu:factor(sex) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome ~ age + edu + factor(sex) + edu:factor(sex)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
- edu:factor(sex) 1 1003253 56582021 3585.3
- age 1 1606902 57185670 3588.4summary(st1) # 最終的に選択されたモデルの結果
Call:
lm(formula = fincome ~ age + edu + factor(sex) + edu:factor(sex),
data = d01)
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -481.771 304.899 -1.580 0.11518
age 9.573 3.312 2.891 0.00414 **
edu 53.000 16.767 3.161 0.00174 **
factor(sex)2 -790.292 375.726 -2.103 0.03630 *
edu:factor(sex)2 61.381 26.874 2.284 0.02310 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
(90 observations deleted due to missingness)
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08names(st1) # オブジェクトに格納されている情報の名前一覧 [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "contrasts" "xlevels" "call"
[13] "terms" "model" "anova" st1$anova # 変数選択の過程で検討されたモデルの分散分析表st1$anovaの結果は，1が出発点のモデル，そこから2，3，4と，AICが低下する限り変数を減少させていっている。
AIC(st1) # ある（単純な）計算方法におけるAIC[1] 4418.357logLik(st1) # 対数尤度，dfとして最尤推定されるパラメタの個数kが示されている。独立変数の個数4個+2になっている。'log Lik.' -2203.179 (df=6)(-2)*logLik(st1)+2*attributes(logLik(st1))$df # -2LL+2kを計算'log Lik.' 4418.357 (df=6)extractAIC(st1) # 別の方法で計算されたAIC[1] 5.000 3582.021step( )関数で使用されているのはextractAIC( )の計算法と同じである事が分かる。AICはモデル同士のAICの差だけが重要で，それについてはextraxtAIC( )で計算してもAIC( )で計算しても同じであるので気にしなくて良い。
# extractAIC(st1)の結果を再現する
log(deviance(st1)/nrow(st1$model))*nrow(st1$model)+2*st1$rank[1] 3582.021変数減少法の選択の最初に与えたモデルについても，AIC( )とextractAIC( )の結果の再現を示しておく。
# step( )の最初のfullモデルについてもextractAIC( )の値を再現しよう
model2 <- lm(fincome ~ age * edu * factor(sex), data = d01)
log(deviance(model2)/nrow(model2$model))*nrow(model2$model)+2*model2$rank[1] 3585.138extractAIC(model2)[1] 8.000 3585.138AIC(model2) # 単純な方のAIC[1] 4421.474(-2)*logLik(model2)+2*attributes(logLik(model2))$df'log Lik.' 4421.474 (df=9)いろいろな変数選択法
さて，漸く，よりカスタマイズした変数選択法の紹介である。
# 変数増加法において，モデルが複雑になるにつれて有効ケースが減少してしまうので，
# 最初から完備ケース分析の準備をしておく必要がある。
full <- complete.cases(d01$fincome, d01$age, d01$edu, d01$sex)
fincome_c <- d01$fincome[full]
age_c <- d01$age[full]
edu_c <- d01$edu[full]
sex_c <- d01$sex[full]
# 完備ケース分析用の変数はデータフレイムに格納していない点に注意
# 年齢＋教育年数＋性別の主効果だけからなるモデルを出発点(object)として，
# それらの二次の交互作用，三次の交互作用まで含むモデルを最大とする(scope)場合の
# 変数増加法(forward)。つまり，交互作用のどれを入れると良いかを検討している
step(object=lm(fincome_c ~ age_c + edu_c + factor(sex_c)),
scope=list(upper=~ age_c * edu_c * factor(sex_c)),
direction="forward")Start: AIC=3585.28
fincome_c ~ age_c + edu_c + factor(sex_c)
Df Sum of Sq RSS AIC
+ edu_c:factor(sex_c) 1 1003253 55578768 3582.0
<none> 56582021 3585.3
+ age_c:edu_c 1 310607 56271414 3585.7
+ age_c:factor(sex_c) 1 16097 56565924 3587.2
Step: AIC=3582.02
fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
+ age_c:edu_c 1 371070 55207698 3582.1
+ age_c:factor(sex_c) 1 1077 55577691 3584.0
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c))
Coefficients:
(Intercept) age_c edu_c
-481.771 9.573 53.000
factor(sex_c)2 edu_c:factor(sex_c)2
-790.292 61.381 # 変数増減法
step(lm(fincome_c ~ age_c * edu_c * factor(sex_c)), direction="both")Start: AIC=3585.14
fincome_c ~ age_c * edu_c * factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c:factor(sex_c) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + age_c:factor(sex_c) +
edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:factor(sex_c) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age_c:edu_c 1 393611 55577691 3584.0
+ age_c:edu_c:factor(sex_c) 1 147674 55036407 3585.1
- edu_c:factor(sex_c) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c 1 371070 55578768 3582.0
<none> 55207698 3582.1
+ age_c:factor(sex_c) 1 23618 55184080 3583.9
- edu_c:factor(sex_c) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
+ age_c:edu_c 1 371070 55207698 3582.1
+ age_c:factor(sex_c) 1 1077 55577691 3584.0
- edu_c:factor(sex_c) 1 1003253 56582021 3585.3
- age_c 1 1606902 57185670 3588.4
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c))
Coefficients:
(Intercept) age_c edu_c
-481.771 9.573 53.000
factor(sex_c)2 edu_c:factor(sex_c)2
-790.292 61.381 # 変数減少法
step(lm(fincome_c ~ age_c * edu_c * factor(sex_c)), direction="backward")Start: AIC=3585.14
fincome_c ~ age_c * edu_c * factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c:factor(sex_c) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + age_c:factor(sex_c) +
edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:factor(sex_c) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age_c:edu_c 1 393611 55577691 3584.0
- edu_c:factor(sex_c) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c 1 371070 55578768 3582.0
<none> 55207698 3582.1
- edu_c:factor(sex_c) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
- edu_c:factor(sex_c) 1 1003253 56582021 3585.3
- age_c 1 1606902 57185670 3588.4
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c))
Coefficients:
(Intercept) age_c edu_c
-481.771 9.573 53.000
factor(sex_c)2 edu_c:factor(sex_c)2
-790.292 61.381 かなり長くなってしまったので，自分でよく見比べて異同を確認して欲しい。変数増加法では，出発点のモデルを変えると，最終結果がbothやbackwardとは異なってしまう事もある。object=lm(fincome_c ~ 1)などとして確認すると良い。
step(object = lm(fincome_c ~ 1),
scope = list(upper = ~ age_c * edu_c * factor(sex_c)),
direction = "forward")Start: AIC=3615.43
fincome_c ~ 1
Df Sum of Sq RSS AIC
+ edu_c 1 5333194 58652205 3591.8
+ age_c 1 867851 63117549 3613.4
<none> 63985400 3615.4
+ factor(sex_c) 1 76310 63909090 3617.1
Step: AIC=3591.85
fincome_c ~ edu_c
Df Sum of Sq RSS AIC
+ age_c 1 1818871 56833335 3584.6
<none> 58652205 3591.8
+ factor(sex_c) 1 283437 58368768 3592.4
Step: AIC=3584.58
fincome_c ~ edu_c + age_c
Df Sum of Sq RSS AIC
<none> 56833335 3584.6
+ age_c:edu_c 1 377115 56456219 3584.6
+ factor(sex_c) 1 251313 56582021 3585.3
Call:
lm(formula = fincome_c ~ edu_c + age_c)
Coefficients:
(Intercept) edu_c age_c
-785.90 75.05 10.16 発展2-2 誤差減少率
前項2-2のst1と同等のモデル（完備ケース分析の処理済みのもの）について，分散説明率と誤差減少率が一致する事を（数式から自明ではあるが）一応確認しておく。
# st1と同じモデル。但し完備ケース分析の処理をしたもの。
reg100203 <- lm(fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c*factor(sex_c))
summary(reg100203) # Multiple R-squaredを確認
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c *
factor(sex_c))
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -481.771 304.899 -1.580 0.11518
age_c 9.573 3.312 2.891 0.00414 **
edu_c 53.000 16.767 3.161 0.00174 **
factor(sex_c)2 -790.292 375.726 -2.103 0.03630 *
edu_c:factor(sex_c)2 61.381 26.874 2.284 0.02310 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08SST <- var(fincome_c)*(length(fincome_c)-1); SST # 従属変数の偏差平方和[1] 63985400SSr <- var(reg100203$residuals)*(length(reg100203$residuals)-1); SSr # 残差平方和[1] 55578768# 誤差減少率の定義式に代入
1-SSr/SST # 平方和から計算[1] 0.13138361-var(reg100203$residuals)/var(fincome_c) # （不偏）分散から計算[1] 0.1313836第10章の【練習問題】の解答
1) の答え
以下の通りで，vifは非常に大きいが，ageとage^2はいずれも5%有意である。
d01$male <- c(1, 0)[d01$sex]
summary(reg10a1 <- lm(edu ~ male + age + I(age^2), data = d01))
Call:
lm(formula = edu ~ male + age + I(age^2), data = d01)
Residuals:
Min 1Q Median 3Q Max
-5.6913 -1.4822 0.2537 1.7052 4.2831
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.571921 3.322319 6.794 4.34e-11 ***
male 0.318687 0.202483 1.574 0.1164
age -0.343527 0.148608 -2.312 0.0213 *
I(age^2) 0.003216 0.001628 1.975 0.0490 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.935 on 372 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.06148, Adjusted R-squared: 0.05392
F-statistic: 8.124 on 3 and 372 DF, p-value: 2.974e-05car::vif(reg10a1) # carパッケイジから関数vif( )だけを呼び出すやり方 male age I(age^2)
1.008407 140.787202 140.688687 age^2の係数が正なので，年齢と教育年数の関係は下に凸な放物線となる。その頂点のx座標（年齢）は，以下の計算から約53.4歳となる。データに含まれるのは30歳～59歳であるが，53歳以下の対象者では，若いほど教育年数が長い傾向があるが，50代半ば以降では（何故か）年長者の方が教育年数が長いという関係になる（50代後半では高学歴層の回答率が高かった可能性も考えられる）。
coef_a1 <- summary(reg10a1)$coefficients
coef_a1 Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.571920597 3.322318546 6.794027 4.340435e-11
male 0.318686634 0.202482838 1.573895 1.163619e-01
age -0.343527205 0.148608268 -2.311629 2.134385e-02
I(age^2) 0.003215804 0.001628121 1.975163 4.898858e-02coef_a1[3, 1]/((-2)*coef_a1[4, 1])[1] 53.412332) の答え
まず問題文の前半の結果は以下の通りである。vifの結果はどれもかなり悪いが，検定は5%水準～10%水準で有意になっている。
summary(reg10a2 <- lm(edu ~ male + age + I(age^2) + male:age, data = d01))
Call:
lm(formula = edu ~ male + age + I(age^2) + male:age, data = d01)
Residuals:
Min 1Q Median 3Q Max
-5.9341 -1.5716 0.3187 1.6643 4.4892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.079177 3.395222 6.208 1.43e-09 ***
male 2.587744 1.169656 2.212 0.0275 *
age -0.300924 0.149609 -2.011 0.0450 *
I(age^2) 0.002997 0.001626 1.844 0.0660 .
male:age -0.049555 0.025162 -1.969 0.0496 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.927 on 371 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.0712, Adjusted R-squared: 0.06118
F-statistic: 7.11 on 4 and 371 DF, p-value: 1.594e-05car::vif(reg10a2)there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif male age I(age^2) male:age
33.90969 143.79339 141.34738 33.95809 問題の後半である。表示の仕方に少し工夫をしたものも付加してある。自分で調べて解読すると良い。
names(reg10a2) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "xlevels" "call" "terms"
[13] "model" coef10a2 <- reg10a2$coefficients # 定数項と偏回帰係数の行列を取り出す
p.age <- 50 # 予測式で求める年齢を，調査対象の30～59歳から自由に指定。
p.male <- 1 # 男性なら1，女性なら0とする
profile <- c(1, p.male, p.age, p.age^2, p.male*p.age) # 独立変数行列
p.male; p.age[1] 1[1] 50sum(coef10a2*profile)[1] 13.63609# 一目でわかる出力を工夫してみる
sprintf("%s，%d歳の教育年数の予測値%.1f年", c("女性", "男性")[p.male+1], p.age, sum(coef10a2*profile))[1] "男性，50歳の教育年数の予測値13.6年"3) の答え
まずは分析結果を確認し，予測式を男女別に数式表現する。
summary(lm(income ~ age + sex + edu, data = d01))
Call:
lm(formula = income ~ age + sex + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 81.298 163.342 0.498 0.61901
age 5.366 1.929 2.782 0.00571 **
sex -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16 予測式は，\(income = 81.3 + 5.4*age - 353*sex + 42*edu\) となる。
ここで，age年齢は調査時満年齢，eduは教育年数であるが，sexは男性=1，女性=2という値である事に注意する。男性と女性でそれぞれ予測式を書き分けると以下の様になる。
男性の予測式：\(income = 81.3 + 5.4*age - 353*1 + 42*edu = -272 + 5.4*age + 42*edu\)
女性の予測式：\(income = 81.3 + 5.4*age - 353*2 + 42*edu = -625 + 5.4*age + 42*edu\)
summary(lm(income ~ age + factor(sex) + edu, data = d01))
Call:
lm(formula = income ~ age + factor(sex) + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
factor(sex)2 -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16factor(sex)2というのは，factor(sex)が2の値を取る場合に，という意味であり，factor(sex)が1の値（男性）の場合は関係しない。
予測式は，\(income = -272 + 5.4*age - 353*factor(sex)2 + 42*edu\) となる。
factor(sex)2の項の意味に注意して予測式を書き分ける。
男性の予測式：\(income = -272 + 5.4*age + 42*edu\)
女性の予測式：\(income = -272 + 5.4*age - 353 + 42*edu = -625 + 5.4*age + 42*edu\)
となり，上の結果と一致する。factor(sex)を投入するのは，女性ダミー変数を投入するのと全く同じ事になる。
d01$female <- c(0, 1)[d01$sex]
summary(lm(income ~ age + female + edu, data = d01))
Call:
lm(formula = income ~ age + female + edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
female -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16定数項に意味を持たせるには
因みに，年齢や教育年数を平均からの偏差に変換すると，定数項は平均的なプロファイルの人の年収平均を意味するので分かりやすい。
d01$d_age <- d01$age - mean(d01$age, na.rm = T)
d01$d_edu <- d01$edu - mean(d01$edu, na.rm = T)
mean(d01$age, na.rm = T); mean(d01$edu, na.rm = T)[1] 45.80208[1] 13.92021summary(lm(income ~ d_age + female + d_edu, data = d01))
Call:
lm(formula = income ~ d_age + female + d_edu, data = d01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 559.267 23.430 23.870 < 2e-16 ***
d_age 5.366 1.929 2.782 0.00571 **
female -352.892 30.670 -11.506 < 2e-16 ***
d_edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16Footnotes
但し独立変数間の相関は極めて強い↩︎
係数の観点からは線形・一次式なので、これは非線形回帰とは呼ばない。非線形回帰は推定すべき係数が非線形になっている場合である。↩︎
data.frame名は簡単に d01 とした。↩︎
まずは1次の項と2次の項を同時投入するのが通常である。↩︎
mean以外にもsd(x, na.rm = T)やfunction(x) sum(!is.na(x))など様々な統計量を使える。↩︎
「こんな事も出来る」と云う例示で，普通は余りこんな事はしない方が良い。↩︎
とは言え，eduの標準誤差は大きめか？↩︎
完備ケース分析のスクリプトの書き方を改訂した（2024.12.04）。↩︎