『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第9章　重回帰分析（Ⅰ）

0　全体の章構成

• 第1章　1変量の記述統計の基礎
• 第2章　2変数の関連の記述統計
• 第3章　推測統計の基礎
• 第4章　統計的検定の一般型
• 第5章　2変数の関連の推定と検定
• 第6章　2群の母平均の差のt検定
• 第7章　平均値の差の分散分析(ANOVA)
• 第8章　単回帰分析(SRA)
• 第9章　重回帰分析（Ⅰ）(MRA)
• 第10章　重回帰分析（Ⅱ）(MRA)
• 第11章　主成分分析と因子分析
• 第12章　ロジスティック回帰分析
• 第13章　階層構造のあるデータの分析

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{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{red}{e_{i}}=\color{green}{\hat{y}_{i}}+\color{red}{e_{i}}\;\;\;(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{b_{j}-\color{magenta}{\beta_{j}}}{\hat{\sigma}_{b_{j}}}\;\sim\;t_{(df=n-p-1)}\] ゼロ仮説\(H_{0}:\beta_{j}=0\)を仮定する場合は，下記の様に検定統計量として計算可能になる（0でなくても，定数の値をゼロ仮説として仮定すれば検定統計量は計算可能になる）。 \[\frac{b_{j}}{\hat{\sigma}_{b_{j}}}\;\sim\;t_{(df=n-p-1)}\] また，自由度が十分に大きくてt分布が標準正規分布で近似出来るならば，下記の式が95%の確率で成り立つ。 \[-1.96<\frac{b_{j}-\color{magenta}{\beta_{j}}}{\hat{\sigma}_{b_{j}}}<1.96\] これを変形して，以下の様な，母偏回帰係数の95%信頼区間が求められる。

\[b_{j}-1.96\hat{\sigma}_{b_{j}}<\color{magenta}{\beta_{j}}<b_{j}+1.96\hat{\sigma}_{b_{j}}\]

　本文の，幸福度を，年齢，教育年数，階層帰属意識，世帯年収で説明する重回帰分析を行い(reg0901)，その結果を表示する。データフレイム名はdata01としているので注意せよ。

d01 <- read.csv("practice.csv")

# まずは各変数の度数分布または基本要約統計量を確認
with(d01, table(fincome, useNA="always"))

fincome
   0   25   75  125  200  300  400  500  600  700  800  925 1125 1375 1750 2500 
   5    1    7    8   17   19   28   36   32   27   23   33   31   13   11    6 
<NA> 
  87 

summary(d01$fincome)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0   400.0   600.0   720.1   925.0  2500.0      87 

with(d01, table(q1700, useNA="always"))

q1700
   0    1    2    3    4    5    6    7    8    9   10   11 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    1    5 

summary(d01$q1700)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  0.000   5.000   7.000   6.694   8.000  11.000       5 

with(d01, table(q1502, useNA="always")) # 10段階階層帰属意識

q1502
   1    2    3    4    5    6    7    8    9   10   11 <NA> 
   2    4   21   66  105   68   62   35    6    6    2    7 

with(d01, table(age, useNA="always")) # 調査時点満年齢

age
  30   31   32   33   34   35   36   37   38   39   40   41   42   43   44   45 
   3    8    5    7   14   10   13   19    9   10   14   15    4   17   21   16 
  46   47   48   49   50   51   52   53   54   55   56   57   58   59 <NA> 
  14   13   17   21   13   12   14   16   11   13   13   14   13   15    0 

with(d01, table(edu, useNA="always")) # 教育年数

edu
   9   11   12   14   16   18 <NA> 
  15    3  110  123  113   12    8 

# 必要に応じて欠損値処理（q1700は少しひねった方法）
d01$q1502 <- c(1:10)[d01$q1502]
d01$q1700 <- c(0:10)[d01$q1700 + 1]
table(d01$q1502, useNA="always")

   1    2    3    4    5    6    7    8    9   10 <NA> 
   2    4   21   66  105   68   62   35    6    6    9 

table(d01$q1700, useNA="always")

   0    1    2    3    4    5    6    7    8    9   10 <NA> 
   2    2    4   17   16   65   42   91   81   41   17    6 

# 重回帰分析を，オブジェクト保存を入れ子状にして実行する
sum0901 <- summary(reg0901 <- lm(q1700 ~ age + edu + q1502 + fincome, d01))
sum0901

Call:
lm(formula = q1700 ~ age + edu + q1502 + fincome, data = d01)

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    
q1502       -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
  (94 observations deleted due to missingness)
Multiple R-squared:  0.1688,    Adjusted R-squared:  0.1571 
F-statistic: 14.47 on 4 and 285 DF,  p-value: 9.035e-11

　偏回帰係数は，年齢と教育年数は全く有意ではなく，階層帰属意識(q1502)と世帯収入は1%水準で有意である。
　標本偏回帰係数(Estimate)をその標準誤差(Std. Error)で割った値がt統計量(t value)であり，その両側検定の有意確率が(Pr(>|t|))である。オブジェクトに格納された情報を使って，この事を以下の様に確認してみよう。

names(reg0901) # lm( )関数の結果

 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "na.action"     "xlevels"       "call"          "terms"        
[13] "model"        

names(sum0901) # summary( )関数の結果

 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"    

coef <- sum0901$coefficients; coef # summary( )関数の一部

                 Estimate   Std. Error   t value     Pr(>|t|)
(Intercept)  7.7782403017 1.0847553998  7.170502 6.448518e-12
age         -0.0159390252 0.0130290266 -1.223347 2.222090e-01
edu          0.0668900748 0.0535015026  1.250247 2.122350e-01
q1502       -0.3223251515 0.0630351876 -5.113416 5.814673e-07
fincome      0.0008729785 0.0002286136  3.818576 1.647257e-04

coef[,1]/coef[,2] # 偏回帰係数÷標準誤差＝t統計量

(Intercept)         age         edu       q1502     fincome 
   7.170502   -1.223347    1.250247   -5.113416    3.818576 

reg0901$df.residual # 偏回帰係数のt検定の自由度

[1] 285

# 両側検定の有意確率なので最後に2倍している。
pt(abs(coef[,1]/coef[,2]), df=reg0901$df.residual, lower.tail=F)*2

 (Intercept)          age          edu        q1502      fincome 
6.448518e-12 2.222090e-01 2.122350e-01 5.814673e-07 1.647257e-04 

それぞれの計算結果が何処に対応しているのかよく確認しよう。因みに，この自由度(df=287)のｔ分布の95%区間の限界値はqt(.025, df=reg0901$df.residual, lower.tail=F)もしくはqt(.975, df=reg0901$df.residual)で求めることができ，1.9683226である。

qt(.025, df=reg0901$df.residual, lower.tail=F)

[1] 1.968323

qt(.975, df=reg0901$df.residual)

[1] 1.968323

　95%信頼区間，90%信頼区間をconfint( )関数で求めた後，90%信頼区間の上限と下限をオブジェクトreg0901やsum0901の情報を使って求めてみる。結果がきちんと一致している事を確認して，信頼区間の関数がどの様な計算をしているのかなどについて理解を深めよう。

confint(reg0901, level=.95)

                    2.5 %       97.5 %
(Intercept)  5.6430917295  9.913388874
age         -0.0415843527  0.009706302
edu         -0.0384181421  0.172198292
q1502       -0.4463987361 -0.198251567
fincome      0.0004229932  0.001322964

confint(reg0901, level=.90)

                      5 %         95 %
(Intercept)  5.9881577270  9.568322876
age         -0.0374397556  0.005561705
edu         -0.0213990524  0.155179202
q1502       -0.4263469346 -0.218303368
fincome      0.0004957163  0.001250241

# 90%信頼区間の下限を計算
coef[,1] - qt(.05, df=reg0901$df.residual, lower.tail=F)*coef[,2]

  (Intercept)           age           edu         q1502       fincome 
 5.9881577270 -0.0374397556 -0.0213990524 -0.4263469346  0.0004957163 

# 90%信頼区間の上限を計算
coef[,1] + qt(.05, df=reg0901$df.residual, lower.tail=F)*coef[,2]

 (Intercept)          age          edu        q1502      fincome 
 9.568322876  0.005561705  0.155179202 -0.218303368  0.001250241 

1-2　分散説明率R²の増分の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( )関数を使用する。また，独立変数の平方和の合計を求めるために少し工夫している。

anova(reg0901)

anova(reg0901)$"Sum Sq"

[1]   1.181081  22.929709  99.320171  41.564379 812.387418

p_x <- 4 # 独立変数の数（モデル自由度）を代入する。
SSmodel <- sum(anova(reg0901)$"Sum Sq"[1:p_x]); SSmodel # モデル平方和

[1] 164.9953

SStotal <- sum(anova(reg0901)$"Sum Sq"[1:(p_x + 1)]); SStotal # 総平方和

[1] 977.3828

SSresidual <- sum(anova(reg0901)$"Sum Sq"[p_x + 1])
SSmodel/SStotal # 分散説明率

[1] 0.1688134

sum0901$r.squared # 自動で計算される分散説明率と一致する事を確認

[1] 0.1688134

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 <- sum(anova(reg0901)$Df[1:4]), dfR <- anova(reg0901)$Df[5], dfT <- sum(anova(reg0901)$Df[1:5])) 
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

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)}\] 　Rでの計算の確認は以下の通りである。

sum0901$r.squared # 分散説明率

[1] 0.1688134

sum0901$fstatistic[2] # F統計量の分子自由度

numdf 
    4 

sum0901$fstatistic[3] # F統計量の分母自由度

dendf 
  285 

sum0901$fstatistic[1] # F統計量

   value 
14.47083 

# F統計量とR2の関係式から計算
sum0901$r.squared/sum0901$fstatistic[2]/((1-sum0901$r.squared)/sum0901$fstatistic[3])

   numdf 
14.47083 

部分モデルの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)\] \(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(q1700, age, edu, q1502, fincome))
d05 <- d01[TorF,c("q1700", "age", "edu", "q1502", "fincome")]

sum0902L <- summary(reg0902L <- lm(q1700 ~ age + edu + q1502 + fincome, d05)) # モデルL
sum0902S <- summary(reg0902S <- lm(q1700 ~ age + edu + q1502, d05)) # モデルS
sum0902L

Call:
lm(formula = q1700 ~ age + edu + q1502 + 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    
q1502       -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-11

sum0902S

Call:
lm(formula = q1700 ~ age + edu + q1502, 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 *  
q1502       -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（4.2526205%ポイント）であり，自由度の増加は1である。

# 分散説明率の増分についての，F統計量の計算
F0 <- (sum0902L$r.squared - sum0902S$r.squared)/((1-sum0902L$r.squared)/sum0902L$df[2]) # F統計量
F0

[1] 14.58153

sum0902L$coefficients[5,3] # モデルLにおけるfincomeの偏回帰係数のt統計量

[1] 3.818576

sum0902L$coefficients[5,3]^2 # モデルLにおけるfincomeの偏回帰係数のt統計量の2乗

[1] 14.58153

# F統計量F0に対応する有意確率の計算
pf(F0, 1, 287, lower.tail=F)

[1] 0.0001644983

モデル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 を使用する。

d05$z_q1700 <- (d05$q1700 - mean(d05$q1700))/sd(d05$q1700)
d05$z_age   <- (d05$age - mean(d05$age))/sd(d05$age)
d05$z_edu   <- (d05$edu - mean(d05$edu))/sd(d05$edu)
d05$z_q1502 <- (d05$q1502 - mean(d05$q1502))/sd(d05$q1502)
d05$z_fincome <- (d05$fincome - mean(d05$fincome))/sd(d05$fincome)
sum0903 <- summary(reg0903 <- lm(z_q1700 ~ z_age + z_edu + z_q1502 + z_fincome, d05))
sum0903

Call:
lm(formula = z_q1700 ~ z_age + z_edu + z_q1502 + z_fincome, data = d05)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.68234 -0.60429  0.06439  0.67542  1.89159 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.501e-16  5.391e-02   0.000 1.000000    
z_age       -6.823e-02  5.577e-02  -1.223 0.222209    
z_edu        7.215e-02  5.771e-02   1.250 0.212235    
z_q1502     -2.817e-01  5.508e-02  -5.113 5.81e-07 ***
z_fincome    2.225e-01  5.827e-02   3.819 0.000165 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9181 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-11

reg0903$coefficients[5] # 世帯年収fincomeの標準化偏回帰係数

z_fincome 
0.2225274 

sd(d05$q1700)/sd(d05$fincome)*reg0903$coefficients[5] # 標準化係数と標準偏差から非標準化係数を計算

   z_fincome 
0.0008729785 

reg0902L$coefficients[5] # 世帯年収fincomeの非標準化係数（1-2より）

     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_q1502     -0.39009413 -0.1732459
z_fincome    0.10782350  0.3372314

階層帰属意識q1502と世帯年収fincomeの標準化係数の絶対値は階層帰属意識の方がやや大きかった。信頼区間（の絶対値）を見ても階層帰属意識の方が下限値も上限値も大きいが，しかし信頼区間の絶対値の広がりはかなり重なっており，母集団に於いてはfincomeの効果の方が大きい事もあり得る。

1-4　分析結果の提示方法

　今までの分析結果からまずは，本文の様な，重回帰分析の結果表を作成する為に必要な情報を出力する。1-2からの完備ケース分析の例を使用する。

sum0902L # 重回帰分析（非標準化解）の結果

Call:
lm(formula = q1700 ~ age + edu + q1502 + 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    
q1502       -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-11

round(sum0902L$coefficients[, c(1, 2, 4)], 4) #　数値の小数桁数を指定

            Estimate Std. Error Pr(>|t|)
(Intercept)   7.7782     1.0848   0.0000
age          -0.0159     0.0130   0.2222
edu           0.0669     0.0535   0.2122
q1502        -0.3223     0.0630   0.0000
fincome       0.0009     0.0002   0.0002

round(sum0903$coefficients[, c(1, 4)], 4) # 標準化係数

            Estimate Pr(>|t|)
(Intercept)   0.0000   1.0000
z_age        -0.0682   0.2222
z_edu         0.0722   0.2122
z_q1502      -0.2817   0.0000
z_fincome     0.2225   0.0002

length(d05$q1700) # ケース数

[1] 290

stargazerパッケイジ

　重回帰分析などの分析結果を表にするのに便利なstargazerパッケイジがあるのでここではそれを利用して見る。

# install.packages("stargazer") # 初めてのPCにはinstallが必要。
library(stargazer)

Please cite as: 

 Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

 R package version 5.2.2. https://CRAN.R-project.org/package=stargazer 

stargazer(reg0902L, type="html", no.space=T, single.row=T, digits=4, intercept.bottom=F,
          title="幸福度(q1700)の重回帰分析",  dep.var.labels   = "幸福度(q1700)", align=T,
          covariate.labels=c("（定数）", "年齢", "教育年数", "階層帰属意識", "世帯年収"),
          notes        = "<sup>&lowast;</sup> <I>p</I> < 0.10; <sup>&lowast;&lowast;</sup> <I>p</I> < 0.05; <sup>&lowast;&lowast;&lowast;</sup> <I>p</I> < 0.01", 
          notes.append = FALSE)

幸福度(q1700)の重回帰分析

	Dependent variable:
	幸福度(q1700)
（定数）	7.7782*** (1.0848)
年齢	-0.0159 (0.0130)
教育年数	0.0669 (0.0535)
階層帰属意識	-0.3223*** (0.0630)
世帯年収	0.0009*** (0.0002)
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(reg0902S, reg0902L, type="html", no.space=T, single.row=T, digits=4, intercept.bottom=F,
          title="幸福度(q1700)の重回帰分析",  dep.var.labels   = "幸福度(q1700)", align=T,
          covariate.labels=c("（定数）", "年齢", "教育年数", "階層帰属意識", "世帯年収"), notes.append = FALSE,
          notes        = "<sup>&lowast;</sup> <I>p</I> < 0.10; <sup>&lowast;&lowast;</sup> <I>p</I> < 0.05; <sup>&lowast;&lowast;&lowast;</sup> <I>p</I> < 0.01",
          omit.stat = c("ser", "f"), keep.stat = c("rsq", "adj.rsq", "aic", "n"))

幸福度(q1700)の重回帰分析

	Dependent variable:
	幸福度(q1700)
	(1)	(2)
（定数）	7.3757*** (1.1050)	7.7782*** (1.0848)
年齢	-0.0074 (0.0131)	-0.0159 (0.0130)
教育年数	0.1304** (0.0520)	0.0669 (0.0535)
階層帰属意識	-0.3659*** (0.0634)	-0.3223*** (0.0630)
世帯年収		0.0009*** (0.0002)
Observations	290	290
R2	0.1263	0.1688
Adjusted R2	0.1171	0.1571
Note:	∗ p < 0.10; ∗∗ p < 0.05; ∗∗∗ p < 0.01

（オプションや注記がうまく機能していない部分もある様だ。）

　標準化係数も併記したい場合には表を自作するしかないかも知れない。

memiscパッケイジのmtable( )関数

　他にも類似の関数がある。

library(memisc)

Warning: package 'memisc' was built under R version 4.0.3

Loading required package: lattice

Loading required package: MASS

Attaching package: 'memisc'

The following object is masked from 'package:ggplot2':

    syms

The following objects are masked from 'package:stats':

    contr.sum, contr.treatment, contrasts

The following object is masked from 'package:base':

    as.array

memisc::mtable(reg0902S, reg0902L, summary.stats=c("R-squared", "adj. R-squared", "F", "p", "AIC", "N"))

Calls:
reg0902S: lm(formula = q1700 ~ age + edu + q1502, data = d05)
reg0902L: lm(formula = q1700 ~ age + edu + q1502 + 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)    
  q1502             -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  

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 から引いたものに等しく，モデルによって従属変数の分散がどの程度減少したかを，元の従属変数の分散に対する割合で表した誤差減少率(Proportinal 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-4の重回帰分析の結果表の情報から自由度調整済み分散説明率を計算して確認してみる。

R2 <- .1673
p <- 4
n <- 292
(adjR2 <- 1-(n-1)/(n-p-1)*(1-R2))

[1] 0.1556944

R2 <- .1247
p <- 3
n <- 292
(adjR2 <- 1-(n-1)/(n-p-1)*(1-R2))

[1] 0.1155823

発展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    4

table("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    4

table("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    4

table(d01$D3, d01$d3, useNA="always")

      
         0   1 <NA>
  0    342   0    0
  1      0  38    0
  <NA>   0   0    4

table(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-02

sum0903$coefficients[c(1, 2), 1]　# これを加算すれば非正規雇用の平均になる

(Intercept)          D2 
   804.3103   -172.5536 

sum(sum0903$coefficients[c(1, 2), 1])

[1] 631.7568

sum(sum0903$coefficients[c(1, 3), 1]) # 自営の平均

[1] 630.3571

sum(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.9

summary(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.1866  
4 - 1 == 0  -133.90      76.30  -1.755   0.2145  
---
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.3922
95% family-wise confidence level
 

Linear Hypotheses:
           Estimate  lwr       upr      
2 - 1 == 0 -172.5536 -330.3544  -14.7528
3 - 1 == 0 -173.9532 -401.9596   54.0532
4 - 1 == 0 -133.9022 -316.4202   48.6159

　最後に，独立変数にカテゴリカル変数を指定した線型モデル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.1924

2) の答え

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-16

anova(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.584585

confint(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