Lesson4.Path analysis and more
Kosugitti
2015年11月13日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
野球のデータ，baseball2015.csvをb.dataというオブジェクトに入れているとします。 できてないひとは，次の通りやってみよう。
b.data <- read.csv("baseball2015.csv",na.strings="*")今日使うパッケージ
- lavaan
- semPaths
重回帰分析の復習
回帰分析とは，\[y=a_1 x_1+a_2 x_2+...+b\] という線形関数を当てはめることでした。 たとえば，年齢を身長と体重で予測するとしたら，次のようになります。
reg1 <- lm(age~height+weight,data=b.data)
summary(reg1)##
## Call:
## lm(formula = age ~ height + weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8751 -2.5722 -0.0974 2.7866 14.2368
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 54.30186 10.24212 5.302 3.06e-07 ***
## height -0.18010 0.06869 -2.622 0.00943 **
## weight 0.11466 0.03956 2.898 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.373 on 197 degrees of freedom
## Multiple R-squared: 0.04348, Adjusted R-squared: 0.03377
## F-statistic: 4.478 on 2 and 197 DF, p-value: 0.01254説明変数は名義尺度水準のものでもかまわない
回帰分析の説明変数が名義尺度水準である場合，このモデルは分散分析とおなじことになります。このように，分散分析（やt検定）と重回帰分析をまとめて一般に一般線形モデルと言います。
確認してみましょう。従属変数は身長だとして，投げる手(throw)を説明変数にすると，次のようになります。
str(b.data$throw) # throw変数はfactor型になっていることに注意！## Factor w/ 2 levels "右","左": 2 1 2 1 1 1 2 1 1 2 ...reg2 <- lm(height~throw,data=b.data)
summary(reg2)##
## Call:
## lm(formula = height ~ throw, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.205 -4.205 -1.205 2.795 22.795
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 182.2047 0.5031 362.183 <2e-16 ***
## throw左 -1.6185 1.3211 -1.225 0.222
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.579 on 198 degrees of freedom
## Multiple R-squared: 0.007523, Adjusted R-squared: 0.00251
## F-statistic: 1.501 on 1 and 198 DF, p-value: 0.222ここで，説明変数に「throw左」とあるのは，throw右(1)にくらべて，左(2)になるとどれぐらい従属変数が上がるか，ということを表しています。そして係数が有意で無い＝傾きが0であるという帰無仮説を棄却できないので，利き手によって身長に差は無い，ということがわかります(これは全く面白くない仮説でしたし！)
これの別表現がt検定ですね。
reg2t <- t.test(height~throw,data=b.data)回帰分析と同じ表現ができていることに注目してください。
同様に，血液型など2水準以上のデータであっても検証できます。
str(b.data$blood) # factor型になっていることを確認## Factor w/ 5 levels "AB型","A型","B型",..: 2 4 2 4 3 4 2 3 2 1 ...reg3 <- lm(height ~ blood,data=b.data)
summary(reg3)##
## Call:
## lm(formula = height ~ blood, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.0667 -3.3875 -0.3833 2.8437 18.6167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 179.8462 1.4967 120.160 < 2e-16 ***
## bloodA型 -0.4628 1.6509 -0.280 0.780
## bloodB型 0.3101 1.7749 0.175 0.861
## bloodO型 0.5538 1.6801 0.330 0.742
## blood不明 9.2205 1.6992 5.426 1.69e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.397 on 195 degrees of freedom
## Multiple R-squared: 0.3423, Adjusted R-squared: 0.3288
## F-statistic: 25.37 on 4 and 195 DF, p-value: < 2.2e-16reg3anova <- aov(height~blood,data=b.data)
summary(reg3anova)## Df Sum Sq Mean Sq F value Pr(>F)
## blood 4 2955 738.7 25.37 <2e-16 ***
## Residuals 195 5679 29.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1回帰分析と同じ数字がでていることを確認しましょう。
パス解析1:予測は続くよどこまでも
重回帰分析は変数Aで変数Bを予測する，という方法ですから，この連鎖をつなげて変数Aが変数Bを，変数Bが変数Cを，変数Cが変数Dを・・・と繋がっていくと「影響力の伝播」を検証できます。 この方法を，パス解析と言います。
例として，年収が体重を，体重が身長を予測するパス解析をやってみましょう。
path1.1 <- lm(weight ~ pay,data=b.data)
path1.2 <- lm(height ~ weight,data=b.data)
summary(path1.1)##
## Call:
## lm(formula = weight ~ pay, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.950 -7.269 -1.624 6.427 35.376
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.517e+01 1.320e+00 64.513 < 2e-16 ***
## pay 2.364e-04 7.843e-05 3.014 0.00291 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.21 on 198 degrees of freedom
## Multiple R-squared: 0.04387, Adjusted R-squared: 0.03905
## F-statistic: 9.086 on 1 and 198 DF, p-value: 0.002913summary(path1.2)##
## Call:
## lm(formula = height ~ weight, data = b.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6964 -2.7573 -0.3232 2.3683 16.1768
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 144.90581 2.49850 58.00 <2e-16 ***
## weight 0.41952 0.02805 14.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.525 on 198 degrees of freedom
## Multiple R-squared: 0.5305, Adjusted R-squared: 0.5281
## F-statistic: 223.7 on 1 and 198 DF, p-value: < 2.2e-16このモデルの背景には，「年収が増えると美味しいものを食べれるから，体重が大きくなるだろう。体重が大きくなったら，足の裏にもお肉がつくから，身長だって伸びるだろう」というストーリー（仮説）があります。
パス解析では，まず仮説がないと分析になりませんし，モデルが当てはまっていても（モデル適合度は\(R^2\)でした)，仮説が実質的に妥当なものでなければ意味がありません。
今回の仮説は，実質的には多少馬鹿馬鹿しいものでしたし，モデルの適合度もよくありませんでしたね。これは悪い例です。みなさんは良い仮説のもと，データに基づいたモデリングを進めてください。
パス解析2:一回で済ませちゃうすごい分析法
パス解析のパスはいくら長くなっても構いませんが，そのたびに回帰分析をやっていくのは面倒です。また，先ほどの例ではPath1.2の適合度はそこそこよかった(\(R^2=0.53\))のですが，Path1.1の適合度はわるかった(\(R^2=0.04\))ですね。この場合，「モデル全体としては，結局よかったの，悪かったの？」ということがはっきりしません。
また，一つのデータに対して，複数のモデルを当てはめることで，一部の変数(体重)が重複して使われることになってしまいました。これはデータの水増しをしているようで，あまり良いことではありません。
こうした問題を解決する手法が構造方程式モデリング(SEM)です。SEMでは，
- 複数のパス経路を一回の分析で行う
- 一回の分析で行うので，モデル適合度も全体に対して算出される
- 潜在変数を設定することができる(次回以降のテーマ)
というメリットがあります。
やってみよう
SEMをRで行う一つの方法は，lavaanパッケージを使うことです。
# install.packages("lavaan") # パッケージが無い場合はインストールをしておく
library(lavaan)## Warning: package 'lavaan' was built under R version 3.2.2## This is lavaan 0.5-19
## lavaan is BETA software! Please report any bugs.# モデルが複数行にわたっていることに注意
model <- '
weight ~ pay
height ~ weight
'
result.path <- sem(model,data=b.data)
summary(result.path,standardized=T,fit.measure=T)## lavaan (0.5-19) converged normally after 17 iterations
##
## Number of observations 200
##
## Estimator ML
## Minimum Function Test Statistic 1.759
## Degrees of freedom 1
## P-value (Chi-square) 0.185
##
## Model test baseline model:
##
## Minimum Function Test Statistic 161.953
## Degrees of freedom 3
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.995
## Tucker-Lewis Index (TLI) 0.986
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3478.834
## Loglikelihood unrestricted model (H1) -3477.955
##
## Number of free parameters 4
## Akaike (AIC) 6965.668
## Bayesian (BIC) 6978.861
## Sample-size adjusted Bayesian (BIC) 6966.189
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.062
## 90 Percent Confidence Interval 0.000 0.210
## P-value RMSEA <= 0.05 0.289
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.026
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Regressions:
## Estimate Std.Err Z-value P(>|z|) Std.lv
## weight ~
## pay 0.000 0.000 3.029 0.002 0.000
## height ~
## weight 0.420 0.028 15.033 0.000 0.420
## Std.all
##
## 0.209
##
## 0.728
##
## Variances:
## Estimate Std.Err Z-value P(>|z|) Std.lv
## weight 124.418 12.442 10.000 0.000 124.418
## height 20.268 2.027 10.000 0.000 20.268
## Std.all
## 0.956
## 0.469回帰分析のときに示された数値に該当するところがどこに出力されているか，見つけておきましょう。ちなみに，回帰分析の推定法(最小二乗法，LS)とちがって，最尤推定法(Maximum Likelihood,LM法)になっているので少し違うことがあります。
さて，出力(summary)の際に二つオプションをつけました。
- standardizedオプションは，標準化係数を出すかどうかを決めます
- fit.measureオプションは適合度指標を出すかどうかを決めます
適合度指標はこのモデルに全体に対して出力されますが，異なる観点から複数の指標が算出されるので，それらを見て総合的に判断します。
| 適合度 | 基準 |
|---|---|
| CFI | 0.9以上なら良いモデル |
| TLI | 0.9以上なら良いモデル |
| AIC | 同じデータセットなら小さいほど良いモデル |
| BIC | 同じデータセットなら小さいほど良いモデル |
| RMSEA | 0.05以下なら良いモデル |
| SRMR | 0.05以下なら良いモデル |
パス解析は，結果を図にして表現するのが一般的です。その方がわかりやすいもんね。 基本的には左から右に，影響力が伝播していくように記載します。 semPlotパッケージを使うと，パス図を自動的に描いてくれるので便利です。
library(semPlot)## Warning: package 'semPlot' was built under R version 3.2.2semPaths(result.path,what="stand",style="lisrel")
本日の課題
- バッターのデータかピッチャーのデータをつかって，重回帰分析の繰り返しによるパス解析を行ってください。
- 同じモデルをlavaanを使って分析してください。
仮説の妥当性，モデルの適合度も評価対象です。良いモデルを作ってくださいね。