d01 <- read.csv("practice.csv")『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第9章 重回帰分析（Ⅰ）
0 全体の章構成
1 重回帰モデル
\[ Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+...+\beta_{j}X_{j}+...+\beta_{p}X_{p}+\epsilon,\;\;\;\;\;\;\epsilon\;\sim\;N(0,\sigma_{\epsilon}^2) \]
\[ \color{blue}y_{i}\color{black}=\color{green}\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\hat{\beta}_{2}x_{2i}+...+\hat{\beta}_{j}x_{ji}+...+\hat{\beta}_{p}x_{pi}\color{black}+\color{red}e_{i}\color{black}=\color{green}\hat{y}_{i}\color{black}+\color{red}e_{i}\color{black}\;\;\;(i=1,2,3,...,n) \]
緑の部分が予測式・予測値である。青が観測値・実測値で，赤が残差である。
1-1 偏回帰係数(partial regression coefficient)のt検定と区間推定
誤差平方の推定値（残差平方和を自由度で割る；残差平均平方）
\[\hat{\sigma}_{\epsilon}^2=\frac{1}{n-p-1}\sum_{i=1}^{n}e_{i}^{2}=\frac{1}{n-p-1}\sum_{i=1}^{n}(y_{i}-\hat{y_{i}})=MS_{residual}\] \(SS_{x_{j}}=\sum_{i=1}^{n}(x_{ji}-\bar{x}_{j})^2\)，
独立変数\(x_{j}\)を他の全ての独立変数\(x_{1},\;x_{2},\;...,x_{p}\)に回帰させた時の分散説明率を\(R_{x_{j}}^2\)とすると，\(b_{j}=\hat{\beta_{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{\color{red}b_{j}\color{black}-\color{blue}\beta_{j}\color{black}}{\color{red}\hat{\sigma}_{b_{j}}\color{black}}\;\sim\;t_{(df=n-p-1)}
\] ゼロ仮説\(H_{0}:\color{blue}\beta_{j}\color{black}=\color{green}0\color{black}\)を仮定すれば，下記の様に検定統計量として計算可能になる。
0でなくても，定数の値をゼロ仮説として仮定すれば検定統計量は計算可能になる。
\[
\frac{\color{red}b_{j}\color{black}}{\color{red}\hat{\sigma}_{b_{j}}\color{black}}\;\sim\;t_{(df=n-p-1)}
\] また，自由度が十分に大きくてt分布が標準正規分布で近似出来るならば，下記の式が95%の確率で成り立つ。 \[
-1.96<\frac{\color{red}b_{j}\color{black}-\color{blue}\beta_{j}\color{black}}{\color{red}\hat{\sigma}_{b_{j}}\color{black}}<1.96
\] これを変形して，以下の様な，母偏回帰係数の95%信頼区間が求められる。
\[\color{red}b_{j}\color{black}-1.96\color{red}\hat{\sigma}_{b_{j}}\color{black}<\color{blue}\beta_{j}\color{black}<\color{red}b_{j}\color{black}+1.96\color{red}\hat{\sigma}_{b_{j}}\color{black}\]
本文の重回帰分析の例
幸福度を，年齢，教育年数，階層帰属意識，世帯年収で説明する重回帰分析を行い(reg0901)，その結果を表示する。データフレイム名はd01。
標本偏回帰係数(Estimate)をその標準誤差(Std. Error)で割った値がt統計量(t value)であり，その両側検定の有意確率が(Pr(>|t|))である。
オブジェクトに格納された情報を使って，これを確認してみよう。
それぞれの計算結果が何処に対応しているのかよく確認しよう。
因みに，この自由度(df=285)のｔ分布の95%区間の限界値はqt(.025, df=reg0901$df.residual, lower.tail = F)もしくはqt(.975, df = reg0901$df.residual)で求めることができる。
qt(.025, df = reg0901$df.residual, lower.tail = F)[1] 1.968323qt(.975, df = reg0901$df.residual)[1] 1.9683231-2 分散説明率R2の増分のF検定
\[ \sum_{i=1}^{n}(y_{i}-\bar{y})^2=\sum_{i=1}^{n}(\hat{y_{i}}-\bar{y}+e_{i})^2=\sum_{i=1}^{n}(\hat{y_{i}}-\bar{\hat{y}})^2+\sum_{i=1}^{n}e_{i}^2 \]
\[ R^2=\frac{SS_{model}}{SS_{total}}=1-\frac{SS_{residual}}{SS_{total}} \]
分散分析表を作成する為に，anova( )関数を使用する。
ano0901 <- anova(reg0901) # type I SSである点に注意。
ano0901分散分析表の作成
ANOVAtable <- matrix(NA, nrow = 3, ncol = 5)
colnames(ANOVAtable) <- c("平方和SS", "自由度df", "平均平方MS", "F値", "p値")
rownames(ANOVAtable) <- c("回帰モデル", "残差", "全体")
ANOVAtable[1:3, 1] <- round(c(SSmodel, SSresidual, SStotal), 2)
ANOVAtable[1:3, 2] <- c(dfM,
dfR <- reg.result$df.residual,
dfT <- dfM + dfR)
ANOVAtable[1:2, 3] <- round(c(SSmodel, SSresidual) / c(dfM, dfR), 3)
Fvalue <- (SSmodel/dfM) / (SSresidual/dfR)
ANOVAtable[1, 5] <- sprintf("%e", pf(Fvalue, dfM, dfR, lower.tail = F))
ANOVAtable[1, 4] <- sprintf("%.3f", Fvalue)
ANOVAtable[3, 3:5] <- ""; ANOVAtable[2, 4:5] <- ""kable(ANOVAtable, align = 'r', caption = "重回帰分析の分散分析表")| 平方和SS | 自由度df | 平均平方MS | F値 | p値 | |
|---|---|---|---|---|---|
| 回帰モデル | 165 | 4 | 41.249 | 14.471 | 9.035426e-11 |
| 残差 | 812.39 | 285 | 2.85 | ||
| 全体 | 977.38 | 289 |
以下で，F値やp値がsum0901の結果と一致している事を確認しよう。
\[ F=\frac{MS_{model}}{MS_{residual}}=\frac{SS_{model}/p}{SS_{residual}/(n-p-1)}\;\sim\;F(p,\;n-p-1) \]
\[ F=\frac{R^2/p}{(1-R^2)/(n-p-1)} \]
部分モデルのF検定
モデルS(Small): \(y=\beta_{0}+\beta_{1}x_{1}+...+\beta_{k}x_{k}+\epsilon\)
モデルL(Large): \(y=\beta_{0}+\beta_{1}x_{1}+...+\beta_{k}x_{k}+\beta_{k+1}x_{k+1}+...+\beta_{p}x_{p}+\epsilon\)
\(H_{0}:\;\beta_{k+1}=\beta_{k+2}=...=\beta_{p}=0\)
\[F=\frac{(R_{L}^2-R_{S}^2)/(p-k)}{(1-R_{L}^2)/(n-p-1)}\sim\;F(p-k, n-p-1)\]
F検定とt検定の同等性
特に \(p=k+1\) の時，
\[F=\frac{R_{L}^2-R_{S}^2}{(1-R_{L}^2)/(n-p-1)}=\left(\frac{b_{p}}{\hat{\sigma}_{b_{p}}}\right)^2=t_{b_{p}}^2\]
偏回帰係数のt検定と（部分）モデルのF検定の同等性について少し詳しく説明しておこう。
Rスクリプトについてはやや異なる工夫を行っている。
# 大きい方のモデルですべての変数が有効であるケースに限定したデータフレイム d05 作成
TorF <- with(d01, complete.cases(happy, age, edu, strata, fincome))
d05 <- d01[TorF, c("happy", "age", "edu", "strata", "fincome")]
sum0902L <- summary(reg0902L <- lm(happy ~ age + edu + strata + fincome, d05)) # モデルL
sum0902S <- summary(reg0902S <- lm(happy ~ age + edu + strata, d05)) # モデルSsum0902L
Call:
lm(formula = happy ~ age + edu + strata + fincome, data = d05)
Residuals:
Min 1Q Median 3Q Max
-4.9329 -1.1113 0.1184 1.2421 3.4787
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.7782403 1.0847554 7.171 6.45e-12 ***
age -0.0159390 0.0130290 -1.223 0.222209
edu 0.0668901 0.0535015 1.250 0.212235
strata -0.3223252 0.0630352 -5.113 5.81e-07 ***
fincome 0.0008730 0.0002286 3.819 0.000165 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.688 on 285 degrees of freedom
Multiple R-squared: 0.1688, Adjusted R-squared: 0.1571
F-statistic: 14.47 on 4 and 285 DF, p-value: 9.035e-11sum0902S
Call:
lm(formula = happy ~ age + edu + strata, data = d05)
Residuals:
Min 1Q Median 3Q Max
-5.2027 -1.2709 0.1112 1.2889 3.7310
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.375711 1.104959 6.675 1.28e-10 ***
age -0.007366 0.013135 -0.561 0.5754
edu 0.130363 0.052047 2.505 0.0128 *
strata -0.365929 0.063447 -5.767 2.09e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.728 on 286 degrees of freedom
Multiple R-squared: 0.1263, Adjusted R-squared: 0.1171
F-statistic: 13.78 on 3 and 286 DF, p-value: 2.031e-08モデルSからモデルLへは分散説明率の増分は0.0425262であり，モデル自由度の増加は1である。
F0 <- (sum0902L$r.squared - sum0902S$r.squared)/((1-sum0902L$r.squared)/sum0902L$df[2]) # F統計量
F0[1] 14.58153# F統計量F0に対応する有意確率の計算
pf(F0, 1, dfR, lower.tail = F)[1] 0.0001648408# モデルLにおけるfincomeの偏回帰係数のt統計量の2乗
sum0902L$coefficients[5,3]^2[1] 14.58153# モデルLにおけるfincomeの偏回帰係数の有意確率
sum0902L$coefficients["fincome", "Pr(>|t|)"][1] 0.0001647257anova(reg0902S, reg0902L)モデルLとモデルSの分散説明率の増加分から計算したF統計量であるF0は，モデルLの独立変数fincomeのt統計量の二乗と一致している。F0のF検定の有意確率もモデルLの独立変数fincomeのt検定の有意確率と一致している。
1-3 標準化偏回帰係数
\(\hat{\beta_{j}}=b_{j}\)，\(\hat{\beta_{j}'}=b_{j}'\)とすると
\[b_{j}=\frac{s_{y}}{s_{x}}b'_{j}\]
1-2の完備ケース分析用データフレイム d05 を使用する。
reg0902L$coefficients["fincome"] # 世帯年収fincomeの非標準化係数（1-2より） fincome
0.0008729785 reg0903$coefficients["z_fincome"] # 世帯年収fincomeの標準化偏回帰係数z_fincome
0.2225274 reg0903$coefficients["z_fincome"] * sd(d05$happy) / sd(d05$fincome) z_fincome
0.0008729785 confint(reg0903, level = .95) 2.5 % 97.5 %
(Intercept) -0.10611408 0.1061141
z_age -0.17800212 0.0415479
z_edu -0.04144032 0.1857443
z_strata -0.39009413 -0.1732459
z_fincome 0.10782350 0.3372314階層帰属意識strataと世帯年収fincomeの標準化係数の絶対値は階層帰属意識の方がやや大きかった。
信頼区間（の絶対値）を見ても階層帰属意識の方が下限値も上限値も大きいが，しかし信頼区間の絶対値の広がりはかなり重なっており，母集団に於いてはfincomeの効果の方が大きい事もあり得る。
1-4 分析結果の提示方法
今までの分析結果からまずは，本文の様な，重回帰分析の結果表を作成する為に必要な情報を出力する。1-2からの完備ケース分析の例を使用する。
knitr::kable(coef.table, align = 'r', caption = "重回帰分析の偏回帰係数の検定")| 偏回帰係数 | s.e. | 標準化係数 | p値 | 95%lower | 95%upper | |
|---|---|---|---|---|---|---|
| (Intercept) | 7.778 | 1.085 | 0.000 | 0.000 | 5.643 | 9.913 |
| z_age | -0.016 | 0.013 | -0.068 | 0.222 | -0.042 | 0.010 |
| z_edu | 0.067 | 0.054 | 0.072 | 0.212 | -0.038 | 0.172 |
| z_strata | -0.322 | 0.063 | -0.282 | 0.000 | -0.446 | -0.198 |
| z_fincome | 0.001 | 0.000 | 0.223 | 0.000 | 0.000 | 0.001 |
カテゴリカル変数も含んだ一般線形モデルの結果表3種
完備ケースデータフレイムで非標準化重回帰
TorF2 <- with(d01, complete.cases(happy, age, EDU2, strata, fincome))
d06 <- d01[TorF2, c("happy", "age", "EDU2", "strata", "fincome")]
sum0902c <- summary(reg0902c <- lm(happy ~ age + EDU2 + strata + fincome, d06))
sum0902c
Call:
lm(formula = happy ~ age + EDU2 + strata + fincome, data = d06)
Residuals:
Min 1Q Median 3Q Max
-5.0018 -1.1529 0.1327 1.2404 3.6865
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.6779115 0.7461300 11.631 < 2e-16 ***
age -0.0183359 0.0129860 -1.412 0.159
EDU2短大卒 0.1995064 0.2443432 0.817 0.415
EDU2四大卒 0.0160403 0.2589339 0.062 0.951
strata -0.3225008 0.0632670 -5.097 6.29e-07 ***
fincome 0.0009644 0.0002289 4.214 3.38e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.693 on 284 degrees of freedom
Multiple R-squared: 0.1667, Adjusted R-squared: 0.152
F-statistic: 11.36 on 5 and 284 DF, p-value: 5.257e-10カテゴリカル変数以外を標準化2
d06$z_happy <- (d06$happy - mean(d06$happy)) /sd(d06$happy)
d06$z_age <- (d06$age - mean(d06$age)) /sd(d06$age)
d06$z_strata <- (d06$strata - mean(d06$strata)) /sd(d06$strata)
d06$z_fincome <- (d06$fincome - mean(d06$fincome))/sd(d06$fincome)sum0903c <- summary(reg0903c <- lm(z_happy ~ z_age + EDU2 + z_strata + z_fincome, d06))
sum0903c
Call:
lm(formula = z_happy ~ z_age + EDU2 + z_strata + z_fincome, data = d06)
Residuals:
Min 1Q Median 3Q Max
-2.71985 -0.62693 0.07216 0.67447 2.00459
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.039398 0.094260 -0.418 0.676
z_age -0.078487 0.055587 -1.412 0.159
EDU2短大卒 0.108486 0.132867 0.817 0.415
EDU2四大卒 0.008722 0.140801 0.062 0.951
z_strata -0.281824 0.055287 -5.097 6.29e-07 ***
z_fincome 0.245820 0.058337 4.214 3.38e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9209 on 284 degrees of freedom
Multiple R-squared: 0.1667, Adjusted R-squared: 0.152
F-statistic: 11.36 on 5 and 284 DF, p-value: 5.257e-10reg.result <- reg0902c # ここに非標準化変数でのlm( )の結果を代入
reg.result.std <- reg0903c # 標準化変数でのlm( )
anova.result <- anova(reg.result)
summary.result <- summary(reg.result)
summary.result.std <- summary(reg.result.std)nVar <- max(reg.result$assign) # 独立変数の数（カテゴリカルは1個と数える）
dfM <- sum(anova.result$Df[1:nVar]) # モデル全体の自由度
SSmodel <- sum(anova.result$"Sum Sq"[1:nVar]) # モデル平方和
SStotal <- sum(anova.result$"Sum Sq"[1:(nVar + 1)]) # 総平方和
SSresidual <- sum(anova.result$"Sum Sq"[nVar + 1])
ANOVAtable <- matrix(NA, nrow = 3, ncol = 5)
colnames(ANOVAtable) <- c("平方和SS", "自由度df", "平均平方MS", "F値", "p値")
rownames(ANOVAtable) <- c("回帰モデル", "残差", "全体")
ANOVAtable[1:3, 1] <- round(c(SSmodel, SSresidual, SStotal), 2)
ANOVAtable[1:3, 2] <- c(dfM,
dfR <- reg.result$df.residual,
dfT <- dfM + dfR)
ANOVAtable[1:2, 3] <- round(c(SSmodel, SSresidual) / c(dfM, dfR), 3)
Fvalue <- (SSmodel/dfM) / (SSresidual/dfR)
ANOVAtable[1, 5] <- sprintf("%e", pf(Fvalue, dfM, dfR, lower.tail = F))
ANOVAtable[1, 4] <- sprintf("%.3f", Fvalue)
ANOVAtable[3, 3:5] <- ""; ANOVAtable[2, 4:5] <- ""PartialRegCoef <- round(summary.result$coefficients, 3) # 数値の小数桁数を指定
StandardizedCoef <- round(summary.result.std$coefficients, 3) # 標準化係数
CI95 <- round(confint(reg.result, level = .95), 3)
coef.table <- matrix(NA, nrow = reg.result$rank, ncol = 6)
colnames(coef.table) <- c("偏回帰係数", "s.e.", "標準化係数", "p値", "95%lower", "95%upper")
rownames(coef.table) <- names(reg.result$coefficients)
coef.table[ , c(1, 2, 4)] <- PartialRegCoef[ , c(1, 2, 4)]
coef.table[ , 3] <- StandardizedCoef[ , 1]
coef.table[ , c(5:6)] <- CI95sup.table <- matrix(NA, nrow = 4, ncol = 1)
rownames(sup.table) <- c("n", "F test p-value ", "R-squared", "adjusted R-squared")
colnames(sup.table) <- ""
sup.table[ , 1] <- c(as.character(nrow(reg.result$model)),
sprintf('%.4e', pf(summary.result$fstatistic[1],
summary.result$fstatistic[2],
summary.result$fstatistic[3], lower.tail = F)),
sprintf('%.4f', summary.result$r.squared),
sprintf('%.4f', summary.result$adj.r.squared))knitr::kable(ANOVAtable, align = 'r', caption = "一般線形モデルの分散分析表")
knitr::kable(coef.table, align = 'r', caption = "一般線形モデルの偏回帰係数の検定")
knitr::kable(sup.table, align = 'r', caption = "分散説明率など")| 平方和SS | 自由度df | 平均平方MS | F値 | p値 | |
|---|---|---|---|---|---|
| 回帰モデル | 162.92 | 5 | 32.583 | 11.362 | 5.256725e-10 |
| 残差 | 814.47 | 284 | 2.868 | ||
| 全体 | 977.38 | 289 |
| 偏回帰係数 | s.e. | 標準化係数 | p値 | 95%lower | 95%upper | |
|---|---|---|---|---|---|---|
| (Intercept) | 8.678 | 0.746 | -0.039 | 0.000 | 7.209 | 10.147 |
| age | -0.018 | 0.013 | -0.078 | 0.159 | -0.044 | 0.007 |
| EDU2短大卒 | 0.200 | 0.244 | 0.108 | 0.415 | -0.281 | 0.680 |
| EDU2四大卒 | 0.016 | 0.259 | 0.009 | 0.951 | -0.494 | 0.526 |
| strata | -0.323 | 0.063 | -0.282 | 0.000 | -0.447 | -0.198 |
| fincome | 0.001 | 0.000 | 0.246 | 0.000 | 0.001 | 0.001 |
| n | 290 |
| F test p-value | 5.2567e-10 |
| R-squared | 0.1667 |
| adjusted R-squared | 0.1520 |
結果をそれぞれcsvファイルに保存する。後で1つのExcelファイルに統合するなどすれば良い。
write.csv(ANOVAtable, file = "MRA_anova.csv", row.names = T)
write.csv(coef.table, file = "MRA_coef.csv", row.names = T)
write.csv(sup.table, file = "MRA_sup.csv", row.names = T)stargazerパッケイジ
重回帰分析などの分析結果を表にするのに便利なstargazerパッケイジがある。
stargazerのコマンド
# install.packages("stargazer") # 初めてのPCにはinstallが必要。
library(stargazer)
stargazer(reg0902L,
type = "html",
no.space = T,
single.row = T,
digits = 4,
title = "幸福度(q1700)の重回帰分析",
dep.var.labels = "幸福度(q1700)",
covariate.labels = c("（定数）", "年齢", "教育年数", "階層帰属意識", "世帯年収"),
notes = "<sup>∗</sup> <I>p</I> < 0.10; <sup>∗∗</sup> <I>p</I> < 0.05; <sup>∗∗∗</sup> <I>p</I> < 0.01",
notes.append = FALSE)| Dependent variable: | |
| 幸福度(q1700) | |
| （定数） | -0.0159 (0.0130) |
| 年齢 | 0.0669 (0.0535) |
| 教育年数 | -0.3223*** (0.0630) |
| 階層帰属意識 | 0.0009*** (0.0002) |
| 世帯年収 | 7.7782*** (1.0848) |
| Observations | 290 |
| R2 | 0.1688 |
| Adjusted R2 | 0.1571 |
| Residual Std. Error | 1.6883 (df = 285) |
| F Statistic | 14.4708*** (df = 4; 285) |
| Note: | ∗ p < 0.10; ∗∗ p < 0.05; ∗∗∗ p < 0.01 |
stargazerでの複数モデル併記コマンド
stargazer(reg0902S, reg0902L,
type = "html",
no.space = T,
single.row = F,
digits = 3,
title = "幸福度(q1700)の重回帰分析",
dep.var.labels = "幸福度(q1700)",
covariate.labels = c("（定数）", "年齢", "教育年数", "階層帰属意識", "世帯年収"),
notes.append = FALSE,
notes = "<sup>∗</sup> <I>p</I> < 0.10; <sup>∗∗</sup> <I>p</I> < 0.05; <sup>∗∗∗</sup> <I>p</I> < 0.01",
omit.stat = c("ser", "f"),
keep.stat = c("rsq", "adj.rsq", "n"))| Dependent variable: | ||
| 幸福度(q1700) | ||
| (1) | (2) | |
| （定数） | -0.007 | -0.016 |
| (0.013) | (0.013) | |
| 年齢 | 0.130** | 0.067 |
| (0.052) | (0.054) | |
| 教育年数 | -0.366*** | -0.322*** |
| (0.063) | (0.063) | |
| 階層帰属意識 | 0.001*** | |
| (0.0002) | ||
| 世帯年収 | 7.376*** | 7.778*** |
| (1.105) | (1.085) | |
| Observations | 290 | 290 |
| R2 | 0.126 | 0.169 |
| Adjusted R2 | 0.117 | 0.157 |
| Note: | ∗ p < 0.10; ∗∗ p < 0.05; ∗∗∗ p < 0.01 | |
複数のモデルの結果を比較する表の一例。
同じデータセットで複数のモデルの優劣を判断する場合，分析に使用されるケースを同一にしておくのが良い。
特に完備ケース分析の場合，それぞれのモデルで有効ケースが異なってしまわないように注意すべきである。
表に整理すべき情報は分析によって異なるかも知れない。非標準化偏回帰係数と標準化係数のいずれを表記すべきか，標準誤差は必要かそれとも有意水準を記すだけで十分かなど，迷う点もあるかも知れない。
memiscパッケイジのmtable( )関数
他にも類似の関数がある。
library(memisc)
memisc::mtable(reg0902S, reg0902L,
summary.stats = c("R-squared", "adj. R-squared", "F", "p", "AIC", "N"))
Calls:
reg0902S: lm(formula = happy ~ age + edu + strata, data = d05)
reg0902L: lm(formula = happy ~ age + edu + strata + fincome, data = d05)
============================================
reg0902S reg0902L
--------------------------------------------
(Intercept) 7.376*** 7.778***
(1.105) (1.085)
age -0.007 -0.016
(0.013) (0.013)
edu 0.130* 0.067
(0.052) (0.054)
strata -0.366*** -0.322***
(0.063) (0.063)
fincome 0.001***
(0.000)
--------------------------------------------
R-squared 0.126 0.169
adj. R-squared 0.117 0.157
F 13.780 14.471
p 0.000 0.000
AIC 1146.183 1133.712
N 290 290
============================================
Significance: *** = p < 0.001;
** = p < 0.01;
* = p < 0.05 偏回帰係数の可視化
ggplot用にdata.frame化
coef0902L <- data.frame(sum0902L$coefficients[, 1], confint(reg0902L))
names(coef0902L) <- c("coefficient", "lower95", "upper95")
coef0902L$Factor <- rownames(sum0902L$coefficients)geom_pointrange( )で信頼区間描画
ggplot(coef0902L[2:nrow(coef0902L), ],
aes(x = Factor, coefficient,
y = coefficient, ymin = lower95, ymax = upper95)) +
geom_hline(yintercept = 0, color = "gray40", linetype = "dashed") +
geom_pointrange(linewidth = 1.2, alpha = .7, shape = 4) +
coord_flip() +
labs(x = NULL,
y = "95% CI of partial regression coefficients",
title = "幸福度の重回帰分析の偏回帰係数の図示",
caption = "（fincomeは単位が「万円」の為，係数の絶対値が非常に小さくなってしまっている。）") +
theme_classic()
1-5 自由度調整済み分散説明率(Adjusted R-squared)
通常の分散説明率（決定係数） \[ R^2=1-\frac{SS_{residual}}{SS_{total}}=1-\frac{\frac{1}{n}SS_{residual}}{\frac{1}{n}SS_{total}}=\frac{\frac{1}{n}SS_{total}-\frac{1}{n}SS_{residual}}{\frac{1}{n}SS_{total}} \] 残差と従属変数の（標本における）分散比を 1 から引いたものに等しく，モデルによって従属変数の分散がどの程度減少したかを，元の従属変数の分散に対する割合で表した誤差減少率(Proportional Reduction in Error; PRE)である。
しかしこれらの分散は標本における統計量であり，母集団に於けるそれらの不偏推定量ではない。つまり母集団に於ける誤差減少率の推定量にはならない。
母集団に於ける誤差減少率の推定量とするには，従属変数の分散を不偏分散に，残差分散を残差の平均平方に書き換えなければならない。そうして得られた統計量が「自由度調整済み分散説明率」である。
\[ adjusted\;R^2=1-\frac{\frac{1}{n-p-1}SS_{residual}}{\frac{1}{n-1}SS_{total}} \]
自由度調整済み分散説明率と無調整の分散説明率の関係は， \[
\frac{SS_{residual}}{SS_{total}}=\frac{n-p-1}{n-1}(1-adjusted\;R^2)=1-R^2
\] より， \[
adjusted\;R^2=1-\frac{n-1}{n-p-1}(1-R^2)
\] である。
自由度調整済み分散説明率は，定義上，無調整の分散説明率よりも値が小さくなる（自分で考えてみよう）。
また，無調整の分散説明率は決して負の値を取らないが，自由度調整済み分散説明率は，よほど当てはまりが悪くて独立変数の数が多い重回帰分析では，負の値を取る事もある。
例えば，n=200，p=10 の時，無調整の決定係数が.05 以下だと自由度調整済み分散説明率はマイナスになる。
これまでの分析結果例でも，自由度調整済み分散説明率が無調整説明率よりも必ず小さな値になっている事が確認出来るだろう。
発展1 1要因分散分析とダミー変数を用いた重回帰分析
第7章発展1では，2群の母平均の差のt検定と1要因分散分析が同値である事を，第8章発展1では，2群の母平均の差のt検定と単回帰分析が同一である事を示した。ここから，要因（独立変数）が2カテゴリの場合には，1要因分散分析と単回帰モデルも同値である事が分かる。本章の発展では，要因（独立変数）が3カテゴリ以上の場合に，一要因分散分析と（複数のダミー変数を用いた）重回帰分析が，同じ一般線型モデル（一般線形モデル）として同値である事を分析例で確認する。
まずはダミー変数の作成方法を紹介しよう。本文では省略したが，最も初歩的な方法だと（例えば）次の様に作る事が出来る。これは殆ど何の工夫もしていない。
# ダミー変数の初期化。代入するものはjobでなくても構わない。
d01$d2 <- d01$d3 <- d01$d4 <- d01$job
d01$d2[d01$job == 1] <- 0; d01$d3[d01$job == 1] <- 0; d01$d4[d01$job == 1] <- 0
d01$d2[d01$job == 2] <- 1; d01$d3[d01$job == 2] <- 0; d01$d4[d01$job == 2] <- 0
d01$d2[d01$job == 3] <- 0; d01$d3[d01$job == 3] <- 1; d01$d4[d01$job == 3] <- 0
d01$d2[d01$job == 4] <- 0; d01$d3[d01$job == 4] <- 0; d01$d4[d01$job == 4] <- 1
table("job"=d01$job, d01$d2, useNA="always") # 意図通りに正しくダミー変数が作成出来たか確認
job 0 1 <NA>
1 182 0 0
2 0 92 0
3 38 0 0
4 68 0 0
<NA> 0 0 4table("job"=d01$job, d01$d3, useNA="always")
job 0 1 <NA>
1 182 0 0
2 92 0 0
3 0 38 0
4 68 0 0
<NA> 0 0 4table("job"=d01$job, d01$d4, useNA="always")
job 0 1 <NA>
1 182 0 0
2 92 0 0
3 38 0 0
4 0 68 0
<NA> 0 0 4本文で紹介したものはこれより少し工夫を施している。上と比べる為にそれぞれD2, D3, D4としよう。
x1 <- c(0, 1, 0, 0) # 非正規ダミー用ヴェクトル
x2 <- c(0, 0, 1, 0) # 自営ダミー用ヴェクトル
x3 <- c(0, 0, 0, 1) # 無職ダミー用ヴェクトル
d01$D2 <- x1[d01$job]; d01$D3 <- x2[d01$job]; d01$D4 <- x3[d01$job]
table(d01$D2, d01$d2, useNA="always") # 上のダミー変数と同一かどうか確認
0 1 <NA>
0 288 0 0
1 0 92 0
<NA> 0 0 4table(d01$D3, d01$d3, useNA="always")
0 1 <NA>
0 342 0 0
1 0 38 0
<NA> 0 0 4table(d01$D4, d01$d4, useNA="always")
0 1 <NA>
0 312 0 0
1 0 68 0
<NA> 0 0 4\[y_{i}=\beta_{0}+\beta_{1}D_{2i}+\beta_{2}D_{3i}+\beta_{3}D_{4i}+e_{i}\]
第7章の分散分析の結果を再確認しておく。
# 第7章の分散分析の結果の確認
tapply(d01$fincome, d01$job, function(x) sum(!is.na(x))) 1 2 3 4
145 74 28 49 tapply(d01$fincome, d01$job, summary)$`1`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
75.0 500.0 700.0 804.3 925.0 2500.0 37
$`2`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 325.0 550.0 631.8 925.0 1750.0 18
$`3`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
75.0 375.0 500.0 630.4 800.0 2500.0 10
$`4`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 200.0 600.0 670.4 925.0 2500.0 19 anova(aov(d01$fincome ~ factor(d01$job)))これを，ダミー変数重回帰で再現してみよう。後の事を考えて，分析結果を一つ一つオブジェクトに格納しながら進めている。
# ダミー重回帰
sum0903 <- summary(reg0903 <- lm(fincome ~ D2 + D3 + D4, d01))
sum0903
Call:
lm(formula = fincome ~ D2 + D3 + D4, data = d01)
Residuals:
Min 1Q Median 3Q Max
-729.3 -304.3 -104.3 190.9 1869.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 804.31 38.35 20.976 < 2e-16 ***
D2 -172.55 65.97 -2.616 0.00936 **
D3 -173.95 95.31 -1.825 0.06901 .
D4 -133.90 76.30 -1.755 0.08031 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 461.7 on 292 degrees of freedom
(88 observations deleted due to missingness)
Multiple R-squared: 0.03038, Adjusted R-squared: 0.02042
F-statistic: 3.05 on 3 and 292 DF, p-value: 0.02898 最下行のF検定の結果は分散分析の結果と一致している。定数項は804.31（万円）となっているが，これが参照カテゴリである正規雇用者の世帯年収の平均と一致し，非正規雇用ダミーD2の偏回帰係数は，非正規雇用の平均631.8と正規雇用の平均804.3との差を示している。
全てのカテゴリの平均を，分析結果を利用して計算・表示しよう。
# 各カテゴリの平均を計算
names(sum0903) [1] "call" "terms" "residuals" "coefficients"
[5] "aliased" "sigma" "df" "r.squared"
[9] "adj.r.squared" "fstatistic" "cov.unscaled" "na.action" sum0903$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 804.3103 38.34513 20.975552 3.220976e-60
D2 -172.5536 65.96544 -2.615818 9.364189e-03
D3 -173.9532 95.31347 -1.825064 6.901289e-02
D4 -133.9022 76.29799 -1.754990 8.030980e-02sum0903$coefficients[c(1, 2), 1] # これを加算すれば非正規雇用の平均になる(Intercept) D2
804.3103 -172.5536 sum(sum0903$coefficients[c(1, 2), 1])[1] 631.7568sum(sum0903$coefficients[c(1, 3), 1]) # 自営の平均[1] 630.3571sum(sum0903$coefficients[c(1, 4), 1]) # 無職の平均[1] 670.4082重回帰分析の結果における，各ダミー変数の偏回帰係数のt検定の部分が，有意水準無調整の多重比較（つまり通常の2群でのt検定）に一致する事も以下の様に確認しよう。
pairwise.t.test(d01$fincome, d01$job, p.adj="none")
Pairwise comparisons using t tests with pooled SD
data: d01$fincome and d01$job
1 2 3
2 0.0094 - -
3 0.0690 0.9891 -
4 0.0803 0.6498 0.7145
P value adjustment method: none 有意確率の小数の表示形式が異なるので分かりにくいかも知れないが，この多重比較の表のうち1列目（つまりカテゴリ1と他のカテゴリの比較の検定の有意確率）の数値は，重回帰分析のCoefficientsの”Pr(>|t|)“の列の結果と一致している。つまり，複数のダミー変数の偏回帰係数のt検定は，単純に二群ごとの平均値の比較のt検定を行っている事と同じなのである。異なる点は，通常の多重比較は全てのカテゴリ間の平均値差を検定するのに対して，ダミー変数重回帰では参照カテゴリと他のカテゴリとの平均値差の検定しか行わない点である。例えば上の重回帰例では，非正規と自営の平均値差の有意確率は計算されない。ここからも，どのカテゴリを参照カテゴリとすべきかが重要である事が分かる。
参照カテゴリとの多重比較
ダミー変数重回帰の偏回帰係数のt検定に対応する，有意水準を調整した多重比較の方法に，Dunnett法がある。本文では脚注で簡単に触れただけであるが，ここで実際に実行してみよう。
まずは，multcompパッケイジを新たにインストールするところから始める（出力結果は省略）。
# まずは使用する追加パッケイジをインストール（PCにつき1回）
# install.packages("multcomp", repos="https://cran.ism.ac.jp")
# インストールしたパッケイジを有効化（Rを新たに起動した場合には必要）
library(multcomp)d01$group <- factor(d01$job)
resultD <- glht(aov(fincome ~ group, d01), linfct=mcp(group="Dunnett"))
resultD
General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Linear Hypotheses:
Estimate
2 - 1 == 0 -172.6
3 - 1 == 0 -174.0
4 - 1 == 0 -133.9summary(resultD)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: aov(formula = fincome ~ group, data = d01)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
2 - 1 == 0 -172.55 65.97 -2.616 0.0273 *
3 - 1 == 0 -173.95 95.31 -1.825 0.1865
4 - 1 == 0 -133.90 76.30 -1.755 0.2144
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method) 先頭のカテゴリとそれ以外のカテゴリの比較だけが検定されている事が読み取れる。有意確率は，偏回帰係数のt検定や無調整のカテゴリ間t検定の有意確率よりも3倍近く高くなっている事が分かる。これが，検定の多重性を調整した結果である。
confint( )関数を使って，平均値差の信頼区間の出力も可能である。95%信頼区間が，有意水準5%の仮説検定の結果と対応する事はもう理解出来ただろうか。
confint(resultD, level=.95)
Simultaneous Confidence Intervals
Multiple Comparisons of Means: Dunnett Contrasts
Fit: aov(formula = fincome ~ group, data = d01)
Quantile = 2.3925
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
2 - 1 == 0 -172.5536 -330.3791 -14.7281
3 - 1 == 0 -173.9532 -401.9953 54.0889
4 - 1 == 0 -133.9022 -316.4488 48.6444最後に，独立変数にカテゴリカル変数を指定した線型モデルlm( )の結果がダミー変数重回帰の結果と完全に一致する事を示しておこう。この様にして，同じ事を様々な方法で実行してみる事を通して，それぞれの分析手法についての理解を深めていけると良い。
summary(lm(fincome ~ factor(job), d01))
Call:
lm(formula = fincome ~ factor(job), data = d01)
Residuals:
Min 1Q Median 3Q Max
-729.3 -304.3 -104.3 190.9 1869.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 804.31 38.35 20.976 < 2e-16 ***
factor(job)2 -172.55 65.97 -2.616 0.00936 **
factor(job)3 -173.95 95.31 -1.825 0.06901 .
factor(job)4 -133.90 76.30 -1.755 0.08031 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 461.7 on 292 degrees of freedom
(88 observations deleted due to missingness)
Multiple R-squared: 0.03038, Adjusted R-squared: 0.02042
F-statistic: 3.05 on 3 and 292 DF, p-value: 0.02898第9章の【練習問題】の解答
1) の答え
出力結果は多少省略した。（本文ではデータフレイム名をdataとしているが，ここではd01としているので注意されたい。）
上の本文のスクリプトよりもやや初心者的なスクリプトで相違があるが，同じ事をやるにも具体的には色々なスクリプトの書き方があると云う事を見て取ってほしい。
# 完備ケース分析．変数名を変えて，変換前後の情報を全て残しておく.
# 新変数は，データフレイムd01の中には含まれていない事に注意．
Full <- complete.cases(d01$income, d01$sex, d01$age, d01$edu)
incomeF <- d01$income[Full]
sexF <- d01$sex[Full]
ageF <- d01$age[Full]
eduF <- d01$edu[Full]
# 男性ダミー変数を作成する
male <- c(1, 0)[sexF]
glm01 <- lm(incomeF ~ male + ageF + eduF)
glm01
Call:
lm(formula = incomeF ~ male + ageF + eduF)
Coefficients:
(Intercept) male ageF eduF
-624.486 352.892 5.366 42.031 summary(glm01)
Call:
lm(formula = incomeF ~ male + ageF + eduF)
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 ***
ageF 5.366 1.929 2.782 0.00571 **
eduF 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
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16 モデルのF検定は1%水準で有意で分散説明率は35%弱，全ての独立変数の偏回帰係数も1%有意である。標本偏回帰係数を見ると，女性より男性は350万円余り個人年収が多く，年齢が1歳上だと5.4万円，教育年数が1年長いと42万円個人年収が多い。
独立変数の設定値をヴェクトル化して，予測値を次の様に計算させる事が出来る。
par_reg <- summary(glm01)$coefficients[, 1]
par_reg(Intercept) male ageF eduF
-624.486186 352.892055 5.366085 42.031195 indep <- c(1, 1, 40, 16) # 男性，40歳，四大卒
sum(par_reg*indep) # 個人年収予測値[1] 615.5484# 女性，35歳，四大卒の個人年収予測値
sum(par_reg*c(1, 0, 35, 16))[1] 235.8259# 女性，50歳，高校卒の個人年収予測値
sum(par_reg*c(1, 0, 50, 12))[1] 148.19242) の答え
summary(lm(scale(incomeF) ~ scale(male) + scale(ageF) + scale(eduF)))
Call:
lm(formula = scale(incomeF) ~ scale(male) + scale(ageF) + scale(eduF))
Residuals:
Min 1Q Median 3Q Max
-2.1132 -0.4997 -0.1359 0.3347 5.2657
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.466e-16 4.417e-02 0.000 1.00000
scale(male) 5.122e-01 4.451e-02 11.506 < 2e-16 ***
scale(ageF) 1.262e-01 4.536e-02 2.782 0.00571 **
scale(eduF) 2.476e-01 4.556e-02 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.812 on 334 degrees of freedom
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16モデルの検定結果は1)と全く同じである（そうあるべきである）。標準化偏回帰係数のt検定の結果も，t統計量と有意確率は全く変化しない。係数の値は変化しており，性別＞教育年数＞年齢の順に絶対値が大きい。よって従属変数に対する効果の大きさはこの順に大きいと言える。
3) の答え
summary(lm(scale(incomeF) ~ male + scale(ageF) + scale(eduF)))
Call:
lm(formula = scale(incomeF) ~ male + scale(ageF) + scale(eduF))
Residuals:
Min 1Q Median 3Q Max
-2.1132 -0.4997 -0.1359 0.3347 5.2657
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.42743 0.05771 -7.406 1.07e-12 ***
male 1.03935 0.09033 11.506 < 2e-16 ***
scale(ageF) 0.12617 0.04536 2.782 0.00571 **
scale(eduF) 0.24759 0.04556 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.812 on 334 degrees of freedom
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16 よく見ると，標準化した年齢と教育年数の結果は2)と全く同じである。男性ダミーは標準化していないので，女性よりも男性の方が個人年収が1標準偏差余り高いと云う事が読み取れる。この影響力は，年齢の8.2標準偏差分，教育年数の4.2標準偏差分に相当する。
2)と3)のどちらがより適切かを言うのは難しいが，そもそも0か1の値しかとらない筈のダミー変数を標準化して「男性ダミー変数が1標準偏差大きいと個人年収が0.51標準偏差大きい」と述べても意味するところが明確だとは思えない。異論もあるかも知れないが，ここではダミー変数に関しては標準化しない立場を提案しておこう。
4) の答え
glm01b <- lm(incomeF ~ male + ageF)
summary(glm01b)
Call:
lm(formula = incomeF ~ male + ageF)
Residuals:
Min 1Q Median 3Q Max
-606.2 -168.3 -68.2 122.4 1897.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 53.181 92.799 0.573 0.567
male 369.394 31.792 11.619 <2e-16 ***
ageF 3.112 1.962 1.586 0.114
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 287.2 on 335 degrees of freedom
Multiple R-squared: 0.2887, Adjusted R-squared: 0.2844
F-statistic: 67.98 on 2 and 335 DF, p-value: < 2.2e-16anova(glm01b, glm01)glm01におけるeduの偏回帰係数のt検定の結果を参照する。
summary(glm01)$coefficients[4, 1:4] # eduFの行 Estimate Std. Error t value Pr(>|t|)
4.203119e+01 7.734500e+00 5.434248e+00 1.063014e-07 # anova(glm01b, glm01) のF検定とp値は一致している
# t統計量の2乗を計算する
summary(glm01)$coefficients[4, 3]^2[1] 29.53105# anova( )のF値と一致している5) の答え
男性ダミー変数や従業上の地位というカテゴリカル変数を含む事を考慮し，非標準化解を優先する。参考の為にそれ以外の変数を標準化した解も後で求めておく。
Full <- complete.cases(d01$income, d01$sex, d01$age, d01$edu, d01$job)
incomeF <- d01$income[Full]
maleF <- c(1, 0)[d01$sex[Full]]
ageF <- d01$age[Full]
eduF <- d01$edu[Full]
jobF <- factor(d01$job[Full])
glm01c <- lm(incomeF ~ maleF + ageF + eduF + jobF)
summary(glm01c)
Call:
lm(formula = incomeF ~ maleF + ageF + eduF + jobF)
Residuals:
Min 1Q Median 3Q Max
-452.68 -135.88 -43.08 71.72 1819.48
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -402.680 147.800 -2.724 0.006782 **
maleF 212.922 34.704 6.135 2.43e-09 ***
ageF 6.390 1.802 3.545 0.000448 ***
eduF 35.339 7.365 4.798 2.42e-06 ***
jobF2 -220.057 40.979 -5.370 1.49e-07 ***
jobF3 -65.775 54.876 -1.199 0.231538
jobF4 -300.329 43.079 -6.972 1.71e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 256.6 on 331 degrees of freedom
Multiple R-squared: 0.439, Adjusted R-squared: 0.4288
F-statistic: 43.17 on 6 and 331 DF, p-value: < 2.2e-16 1)の結果と比較すると，分散説明率は35%弱から44%弱に大きく上昇した。男性ダミーの効果は140万円ほど低下している。年齢の効果は微増（10歳で54万円から64万円へ），教育年数の効果は微減（4年で170万円から140万円へ）したが大きな違いではない。新たに追加投入した従業上の地位は正規雇用が参照カテゴリ（基準カテゴリ）になっており，自営は正規より66万円ほど低いが統計的に有意ではない。非正規は正規より220万円有意に低く，無職は正規より300万円有意に低い。男性ダミーの効果の減少分の多くは，従業上の地位の効果に表れたと考えられる。つまり，男女で従業上の地位（特に正規，非正規，無職）の分布が異なり，その影響が最初のモデルでは男性ダミーの係数に反映されていたのだと考えられる。
以下は，ダミー変数以外を標準化した標準化解の結果である。
# 標準化解
glm01d <- lm(scale(incomeF) ~ maleF + scale(ageF) + scale(eduF) + jobF)
summary(glm01d)
Call:
lm(formula = scale(incomeF) ~ maleF + scale(ageF) + scale(eduF) +
jobF)
Residuals:
Min 1Q Median 3Q Max
-1.3333 -0.4002 -0.1269 0.2112 5.3588
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.09137 0.09309 0.982 0.327034
maleF 0.62711 0.10221 6.135 2.43e-09 ***
scale(ageF) 0.15025 0.04238 3.545 0.000448 ***
scale(eduF) 0.20817 0.04338 4.798 2.42e-06 ***
jobF2 -0.64812 0.12069 -5.370 1.49e-07 ***
jobF3 -0.19372 0.16162 -1.199 0.231538
jobF4 -0.88454 0.12688 -6.972 1.71e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7558 on 331 degrees of freedom
Multiple R-squared: 0.439, Adjusted R-squared: 0.4288
F-statistic: 43.17 on 6 and 331 DF, p-value: < 2.2e-16# 係数の信頼区間も求めてみる
confint(glm01c) 2.5 % 97.5 %
(Intercept) -693.426546 -111.933520
maleF 144.653486 281.190717
ageF 2.844591 9.935687
eduF 20.851302 49.826177
jobF2 -300.669729 -139.444016
jobF3 -173.723743 42.174668
jobF4 -385.072714 -215.584585confint(glm01d) 2.5 % 97.5 %
(Intercept) -0.09174970 0.2744974
maleF 0.42603877 0.8281732
scale(ageF) 0.06688499 0.2336183
scale(eduF) 0.12282620 0.2935049
jobF2 -0.88554354 -0.4106956
jobF3 -0.51165755 0.1242144
jobF4 -1.13413032 -0.6349476