例題: 教育年数が職業の社会的地位の高さに与える効果

職業の社会的地位の高さを職業威信スコアによって測定した。教育年数が長いほど、高い社会的地位の高い職に就く傾向があるかを調べたい。ただし年齢の効果を統制して分析したい。

模擬データの作成

 age <- c(62, 70, 58, 50, 58, 23, 64, 31, 43, 68, 65, 21, 38, 30, 28, 60, 33, 57, 57, 32)
 eduy <- c(16, 9, 12, 12, 12, 12, 12, 16, 16, 9, 12, 12, 12, 12, 12, 9, 9, 12, 12, 12)
 job_sc <- c(84.3, 45.6, 52.2, 42.4, 42.0, 48.2, 44.3, 47.2, 52.2, 
             51.3, 39.0, 38.1, 52.2, 52.2, 59.7, 48.9, 39.0, 48.9, 48.9, 47.9)
 dat <- data.frame(
   age = age,             # 年齢 
   eduy = eduy, 　　　　  # 教育年数
   job_sc = job_sc)       # 職業威信スコア
 
# データを見てみる（最初の6行分)
head(dat)
##   age eduy job_sc
## 1  62   16   84.3
## 2  70    9   45.6
## 3  58   12   52.2
## 4  50   12   42.4
## 5  58   12   42.0
## 6  23   12   48.2

分析

  • 従属変数：職業威信スコア（job_sc）
  • 独立変数：教育年数（eduy）、年齢（age）（※年齢は統制変数としての位置づけ）
# モデルに従って計算をし、outというobjectに付ける
out <- lm(formula = job_sc ~ eduy + age,      # 基本的にはobject名で式を書く 
          data    = dat)　　　　　　　　　　　# データ元

# 結果を見てみる
# 教育年数は職業維新スコアを2.4程度上げる効果がある（5%水準で有意）
summary(out)   
## 
## Call:
## lm(formula = job_sc ~ eduy + age, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.705  -7.210  -0.310   4.498  23.427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  13.7739    14.8724   0.926   0.3673  
## eduy          2.3976     1.0068   2.381   0.0292 *
## age           0.1409     0.1292   1.091   0.2906  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.985 on 17 degrees of freedom
## Multiple R-squared:  0.2605, Adjusted R-squared:  0.1735 
## F-statistic: 2.994 on 2 and 17 DF,  p-value: 0.0769

多重共線性の確認

# carパッケージをインストールする（未インストールの時のみ）
# install.packages("car")         # 未インストールの場合のみ#を外してこのコマンドを実行する

# パッケージを読み込む
library(car)
 
# 多重共線性を計算する
vif(out)                # いずれもVIF < 2 で多重共線性の心配はない
##     eduy      age 
## 1.054646 1.054646

誤差分布の確認

plot(out)

1 標準化回帰係数を用いた回帰分析

1.1 方法1: 個別にobjectを作成する

先にscaleを用いて新しいobjectを作成し、それを回帰分析に用いる

# 従属変数と独立変数を標準化する
dat$age.z <- scale(dat$age)
dat$eduy.z <- scale(dat$eduy)
dat$job_sc.z <- scale(dat$job_sc)

out2 <- lm(job_sc.z ~ eduy.z + age.z, data = dat)
summary(out2)

1.2 方法2: 式中にscale()を加える

以下のように式中の変数名を直接scale（）で囲ってもよい。

out3 <- lm(scale(job_sc) ~ scale(eduy) + scale(age), data = dat)
summary(out3)

2 階層的回帰分析と見やすい結果の並べ方

  • model 1: 独立変数は教育年数のみ
  • model 2: 独立変数は教育年数＋年齢（年齢をmodel 1に加える）

update関数の使用

# model 1
md1.out <- lm(formula = job_sc ~ eduy, 
              data    = dat)　　　　　　　　　　　

# model 2:
md2.out <- update(md1.out,          # モデル１の結果オブジェクトを参照する
                  . ~ . + age)　　　# モデル１の式に、ageを加える（モデル1と変わらない部分は.で省略できる）　　　　　　　　

結果を並べる

以上二つの結果を横に並べる

# パッケージのインストール（未インストールの時のみ）
# install.packages("texreg")          

# ライブラリを読み込む
library(texreg)

# スクリーン上で確認する
screenreg(list(md1.out,             # 年齢を統制した場合、教育年数の効果はより強まっていることが確認できる
               md2.out)) 
## 
## ===============================
##              Model 1   Model 2 
## -------------------------------
## (Intercept)   23.45     13.77  
##              (12.00)   (14.87) 
## eduy           2.15 *    2.40 *
##               (0.99)    (1.01) 
## age                      0.14  
##                         (0.13) 
## -------------------------------
## R^2            0.21      0.26  
## Adj. R^2       0.16      0.17  
## Num. obs.     20        20     
## ===============================
## *** p < 0.001; ** p < 0.01; * p < 0.05
# 結果をhtmlで出力することもできる(出力されたファイルは、自分のプロジェクトフォルダに入る)
htmlreg(list(md1.out,
             md2.out),
        custom.model.names = c("1", "2"),                     # モデル名の変更
        custom.coef.names = c("切片", "教育年数", "年齢"),    # 日本語名に変数名を日本語に直す場合
        file = "regression_table.html")